module m_MEF90_Elements
#include "petsc/finclude/petsc.h"
   use m_MEF90_Parameters
   use m_MEF90_Utils
   use m_MEF90_LinAlg
   use petscdmplex
   use, intrinsic :: iso_c_binding
   implicit none(type, external)

   !Private
   public :: MEF90ElementCreate
   public :: MEF90ElementDestroy

   public :: MEF90ElementType
   public :: MEF90Element2DVect, MEF90Element2DScal
   public :: MEF90Element3DVect, MEF90Element3DScal
   public :: MEF90ElementGetType, MEF90ElementGetTypeBoundary

   type MEF90ElementType
      ! name is the element name in english language
      character(len=MEF90MXSTRLEN)                 :: Name

      !!! Number of dof at each
      !!!   vertex edge face cell in 3D,
      !!!   vertex edge cell in 2D (last value is unused, i.e. always use depth to locate DoF)
      !!! i.e. point depth -1
      PetscInt, dimension(4)                        :: numDofs
      PetscInt                                     :: numDoF
      !!! numDof = numVertexDof * numVertex + numEdgeDof * numEdge
      !!!          + numFaceDof * numFace + numCellDof
      PetscInt                                     :: dim
      !!! The dimension of the cell associated with the element
      PetscInt                                     :: coDim
      !!! The co-dimension of the cell associated with the element
      PetscInt                                     :: order
      !! The polynomial order
   end type MEF90ElementType

   enum, bind(c)
      enumerator :: &
         MEF90ElementFamilyLagrange = 0
   end enum

   character(kind=c_char, len=MEF90MXSTRLEN), dimension(4), parameter, public   :: MEF90ElementFamily = [ &
                                                                                   "Lagrange           ", &  ! 0
                                                                                   "MEF90ElementFamily ", &
                                                                                   "prefix_            ", &
                                                                                   c_null_char//"                  "]
   character(len=MEF90MXSTRLEN), dimension(4), protected  :: MEF90ElementFamilyList

   type(MEF90ElementType), parameter, public :: MEF90_NULL_ELEMENT = MEF90ElementType( &
                                                "NULL", &  ! name
                                                [0, 0, 0, 0], 0, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                0, 0, 0 &  ! dim,codim,order
                                                )

   !!!
   !!! P0 pseudo-elements. These are used for creating sections only
   !!!
   type(MEF90ElementType), parameter, public :: MEF90P0Lagrange2D = MEF90ElementType( &
                                                "MEF90P0Lagrange2D", &  ! name
                                                [0, 0, 1, 0], 1, &  ! numVertexDof,numEdgeDof,numCellDof,unused,numDof
                                                2, 0, 0 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90P0Lagrange3D = MEF90ElementType( &
                                                "MEF90P0Lagrange3D", &  ! name
                                                [0, 0, 0, 1], 1, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 0, 0 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90P0Lagrange2DBoundary = MEF90ElementType( &
                                                "MEF90P0Lagrange2DBoundary", &  ! name
                                                [0, 1, 0, 0], 1, &  ! numVertexDof,numEdgeDof,numCellDof,unused,numDof
                                                2, 0, 0 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90P0Lagrange3DBoundary = MEF90ElementType( &
                                                "MEF90P0Lagrange3DBoundary", &  ! name
                                                [0, 0, 1, 0], 1, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 0, 0 &  ! dim,codim,order
                                                )
   !!!
   !!! P1 Simplicial Lagrange elements
   !!!
   type(MEF90ElementType), parameter, public :: MEF90P1Lagrange2D = MEF90ElementType( &
                                                "MEF90P1Lagrange2D", &  ! name
                                                [1, 0, 0, 0], 3, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                2, 0, 1 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90P1Lagrange3D = MEF90ElementType( &
                                                "MEF90P1Lagrange3D", &  ! name
                                                [1, 0, 0, 0], 4, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 0, 1 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90P1Lagrange2DBoundary = MEF90ElementType( &
                                                "MEF90P1Lagrange2DBoundary", &  ! name
                                                [1, 0, 0, 0], 2, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                2, 1, 1 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90P1Lagrange3DBoundary = MEF90ElementType( &
                                                "MEF90P1Lagrange3DBoundary", &  ! name
                                                [1, 0, 0, 0], 3, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 1, 1 &  ! dim,codim,order
                                                )
   !!!
   !!! P2 Simplicial Lagrange elements
   !!!
   type(MEF90ElementType), parameter, public :: MEF90P2Lagrange2D = MEF90ElementType( &
                                                "MEF90P2Lagrange2D", &  ! name
                                                [1, 1, 0, 0], 6, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                2, 0, 2 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90P2Lagrange3D = MEF90ElementType( &
                                                "MEF90P2Lagrange3D", &  ! name
                                                [1, 1, 0, 0], 10, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 0, 2 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90P2Lagrange2DBoundary = MEF90ElementType( &
                                                "MEF90P2Lagrange2DBoundary", &  ! name
                                                [1, 1, 0, 0], 3, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                2, 1, 2 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90P2Lagrange3DBoundary = MEF90ElementType( &
                                                "MEF90P2Lagrange3DBoundary", &  ! name
                                                [1, 1, 0, 0], 6, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 1, 2 &  ! dim,codim,order
                                                )
   !!!
   !!! Q1 tensor product Lagrange elements
   !!!
   type(MEF90ElementType), parameter, public :: MEF90Q1Lagrange2D = MEF90ElementType( &
                                                "MEF90Q1Lagrange2D", &  ! name
                                                [1, 0, 0, 0], 4, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                2, 0, 1 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90Q1Lagrange3D = MEF90ElementType( &
                                                "MEF90Q1Lagrange3D", &  ! name
                                                [1, 0, 0, 0], 8, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 0, 1 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90Q1Lagrange2DBoundary = MEF90ElementType( &
                                                "MEF90Q1Lagrange2DBoundary", &  ! name
                                                [1, 0, 0, 0], 2, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                2, 1, 1 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90Q1Lagrange3DBoundary = MEF90ElementType( &
                                                "MEF90Q1Lagrange3DBoundary", &  ! name
                                                [1, 0, 0, 0], 4, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 1, 1 &  ! dim,codim,order
                                                )
   !!!
   !!! Q2 tensor product Lagrange elements
   !!!
   type(MEF90ElementType), parameter, public :: MEF90Q2Lagrange2D = MEF90ElementType( &
                                                "MEF90Q2Lagrange2D", &  ! name
                                                [1, 1, 0, 1], 9, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                2, 0, 1 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90Q2Lagrange3D = MEF90ElementType( &
                                                "MEF90Q2Lagrange3D", &  ! name
                                                [1, 1, 1, 1], 27, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 0, 1 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90Q2Lagrange2DBoundary = MEF90ElementType( &
                                                "MEF90Q2Lagrange2DBoundary", &  ! name
                                                [1, 1, 0, 0], 3, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                2, 1, 1 &  ! dim,codim,order
                                                )
   type(MEF90ElementType), parameter, public :: MEF90Q2Lagrange3DBoundary = MEF90ElementType( &
                                                "MEF90Q2Lagrange3DBoundary", &  ! name
                                                [1, 1, 0, 1], 9, &  ! numVertexDof,numEdgeDof,numFaceDof,numCellDof,numDof
                                                3, 1, 1 &  ! dim,codim,order
                                                )
   type MEF90Element2DScal
      PetscReal, dimension(:, :), pointer             :: BF => null()
      type(Vect2D), dimension(:, :), pointer          :: Grad_BF => null()
      PetscReal, dimension(:), pointer               :: Gauss_C => null()
      type(Vect2D)                                 :: outerNormal ! only makes sense for elements of codim 1
   end type MEF90Element2DScal

   type MEF90Element2DVect
      type(Vect2D), dimension(:, :), pointer         :: BF => null()
      type(Mat2D), dimension(:, :), pointer          :: Grad_BF => null()
      type(MatS2D), dimension(:, :), pointer         :: GradS_BF => null()
      PetscReal, dimension(:), pointer               :: Gauss_C => null()
      type(Vect2D)                                 :: outerNormal ! only makes sense for elements of codim 1
   end type MEF90Element2DVect

   type MEF90Element3DScal
      PetscReal, dimension(:, :), pointer             :: BF => null()
      type(Vect3D), dimension(:, :), pointer         :: Grad_BF => null()
      PetscReal, dimension(:), pointer               :: Gauss_C => null()
      type(Vect3D)                                 :: outerNormal ! only makes sense for elements of codim 1
   end type MEF90Element3DScal

   type MEF90Element3DVect
      type(Vect3D), dimension(:, :), pointer         :: BF => null()
      type(Mat3D), dimension(:, :), pointer          :: Grad_BF => null()
      type(MatS3D), dimension(:, :), pointer         :: GradS_BF => null()
      PetscReal, dimension(:), pointer               :: Gauss_C => null()
      type(Vect3D)                                 :: outerNormal ! only makes sense for elements of codim 1
   end type MEF90Element3DVect

   interface MEF90ElementCreate
      module procedure Element2DScalInitSet, Element2DVectInitSet, &
         Element3DScalInitSet, Element3DVectInitSet
   end interface MEF90ElementCreate

   interface MEF90ElementDestroy
      module procedure Element2DScalDestroy, Element2DVectDestroy, &
         Element3DScalDestroy, Element3DVectDestroy, &
         Element2DScalDestroySet, Element2DVectDestroySet, &
         Element3DScalDestroySet, Element3DVectDestroySet
   end interface MEF90ElementDestroy

contains
#undef __FUNCT__
#define __FUNCT__ "MEF90ElementsInitialize_Private"
!!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca)
!!!
!!!  MEF90ElementsInitialize_Private:
!!!

   subroutine MEF90ElementsInitialize_Private(ierr)
      PetscErrorCode, intent(OUT)                   :: ierr

      MEF90ElementFamilyList(1) = 'Lagrange'
      MEF90ElementFamilyList(2) = 'MEF90ElementFamily'
      MEF90ElementFamilyList(3) = '_MEF90ElementFamily'
      MEF90ElementFamilyList(4) = ''
      ierr = 0
   end subroutine MEF90ElementsInitialize_Private

#undef __FUNCT__
#define __FUNCT__ "MEF90ElementGetType"
!!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca)
!!!
!!!  MEF90ElementGetType: Return an element given a family and order
!!!

   subroutine MEF90ElementGetType(elemFamily, order, cellType, elemType, ierr)
      PetscEnum, intent(IN)                             :: elemFamily
      PetscInt, intent(IN)                              :: order
      type(eDMPolytopeType), intent(IN)                 :: cellType
      type(MEF90ElementType), intent(OUT)               :: elemType
      PetscErrorCode, intent(INOUT)                     :: ierr

      elemType = MEF90_NULL_ELEMENT
      select case (elemFamily)
      case (MEF90ElementFamilyLagrange)
         select PetscEnumCase(cellType)
         PetscEnumCase(DM_POLYTOPE_TRIANGLE)
         select case (order)
         case (0)
            elemType = MEF90P0Lagrange2D
         case (1)
            elemType = MEF90P1Lagrange2D
         case (2)
            elemType = MEF90P2Lagrange2D
            ! Case Default
            !    Write(*,*) __FUNCT__,': Unimplemented order',order
            !    ierr = PETSC_ERR_SUP
         end select ! order
         PetscEnumCase(DM_POLYTOPE_TETRAHEDRON)
         select case (order)
         case (0)
            elemType = MEF90P0Lagrange3D
         case (1)
            elemType = MEF90P1Lagrange3D
         case (2)
            elemType = MEF90P2Lagrange3D
            ! Case Default
            !    Write(*,*) __FUNCT__,': Unimplemented order',order
            !    ierr = PETSC_ERR_SUP
         end select ! order
         PetscEnumCase(DM_POLYTOPE_QUADRILATERAL)
         select case (order)
         case (0)
            elemType = MEF90P0Lagrange2D
         case (1)
            elemType = MEF90Q1Lagrange2D
         case (2)
            elemType = MEF90Q2Lagrange2D
            ! Case Default
            !    Write(*,*) __FUNCT__,': Unimplemented order',order
            !    ierr = PETSC_ERR_SUP
         end select ! order
         PetscEnumCase(DM_POLYTOPE_HEXAHEDRON)
         select case (order)
         case (0)
            elemType = MEF90P0Lagrange3D
         case (1)
            elemType = MEF90Q1Lagrange3D
         case (2)
            elemType = MEF90Q2Lagrange3D
            ! Case Default
            !    Write(*,*) __FUNCT__,': Unimplemented order',order
            !    ierr = PETSC_ERR_SUP
         end select ! order
         ! Case Default
         !    Write(*,*) __FUNCT__,': Unknown cell type (see $PETSC_DIR/srcdm/f90-mod/petscdm.h)',cellType
         !    ierr = PETSC_ERR_SUP
      end select ! cellType
      ! Case Default
      ! Write(*,*) __FUNCT__,': Unknown element family',elemFamily
      ! ierr = PETSC_ERR_SUP
      end select
      end subroutine MEF90ElementGetType

#undef __FUNCT__
#define __FUNCT__ "MEF90ElementGetTypeBoundary"
!!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca)
!!!
!!!  MEF90ElementGetTypeBoundary: Return an element given a family and order
!!!

      subroutine MEF90ElementGetTypeBoundary(elemFamily, order, cellType, elemType, ierr)
         PetscEnum, intent(IN)                             :: elemFamily
         PetscInt, intent(IN)                              :: order
         type(eDMPolytopeType), intent(IN)                 :: cellType
         type(MEF90ElementType), intent(OUT)               :: elemType
         PetscErrorCode, intent(INOUT)                     :: ierr

         elemType = MEF90_NULL_ELEMENT
         select case (elemFamily)
         case (MEF90ElementFamilyLagrange)
            select PetscEnumCase(cellType)
            PetscEnumCase(DM_POLYTOPE_SEGMENT)
            select case (order)
            case (0)
               elemType = MEF90P0Lagrange2DBoundary
            case (1)
               elemType = MEF90P1Lagrange2DBoundary
            case (2)
               elemType = MEF90P2Lagrange2DBoundary
               ! Case Default
               !    Write(*,*) __FUNCT__,': Unimplemented order',order
               !    ierr = PETSC_ERR_SUP
            end select ! order
            PetscEnumCase(DM_POLYTOPE_TRIANGLE)
            select case (order)
            case (0)
               elemType = MEF90P0Lagrange3DBoundary
            case (1)
               elemType = MEF90P1Lagrange3DBoundary
            case (2)
               elemType = MEF90P2Lagrange3DBoundary
               ! Case Default
               !    Write(*,*) __FUNCT__,': Unimplemented order',order
               !    ierr = PETSC_ERR_SUP
            end select ! order
            PetscEnumCase(DM_POLYTOPE_QUADRILATERAL)
            select case (order)
            case (0)
               elemType = MEF90P0Lagrange3DBoundary
            case (1)
               elemType = MEF90Q1Lagrange3DBoundary
            case (2)
               elemType = MEF90Q2Lagrange3DBoundary
               ! Case Default
               !    Write(*,*) __FUNCT__,': Unimplemented order',order
               !    ierr = PETSC_ERR_SUP
            end select ! order
            ! Case Default
            !    Write(*,*) __FUNCT__,': Unknown cell type (see $PETSC_DIR/src/dm/f90-mod/petscdm.h)',cellType
            !    ierr = PETSC_ERR_SUP
         end select ! cellType
         ! Case Default
         ! Write(*,*) __FUNCT__,': Unknown element family',elemFamily
         ! ierr = PETSC_ERR_SUP
         end select
      end subroutine MEF90ElementGetTypeBoundary

#undef __FUNCT__
#define __FUNCT__ "Element2DScalInitSet"
!!! author: Blaise Bourdin (2014-2022, bourdin@lsu.edu)
!!! author: Blaise Bourdin (2022, bourdin@mcmaster.ca)
!!!
!!!  Element2DScalInitSet:
!!!

      subroutine Element2DScalInitSet(dm, cellIS, dElem, dQuadratureOrder, elemType, ierr)
         type(tDM), intent(IN)                             :: dm
         type(tIS), intent(IN)                             :: cellIS
         type(MEF90Element2DScal), dimension(:), pointer    :: dElem
         PetscInt, intent(IN)                              :: dQuadratureOrder
         type(MEF90ElementType), intent(IN)                :: elemType
         PetscErrorCode, intent(OUT)                       :: ierr

         PetscInt                                         :: iELoc
         PetscInt, dimension(:), pointer                    :: CellID
         PetscReal, dimension(:), pointer                 :: v0, BB, BBinv
         PetscReal                                        :: detBBinv
         type(Mat2D)                                      :: Bt
         PetscReal                                        :: length
         PetscReal, dimension(:), pointer                   :: centroid, innerNormal
         type(Vect2D)                                     :: outerNormal

         select case (elemType%name)
         case (MEF90P1Lagrange2D%name, MEF90P2Lagrange2D%name)
            PetscCall(ISGetIndices(CellIS, CellID, ierr))
            allocate (dElem(size(cellID)), stat=ierr)
            if (size(CellID) > 0) then
               allocate (v0(2))
               allocate (BB(4))
               allocate (BBinv(4))
               do iELoc = 1, size(CellID)
                  PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cellID(iELoc), v0, BB, BBinv, detBBinv, ierr))
                  !! Petsc uses a reference simplex with vertices at (-1,-1), (1,-1) and (-1,1)
                  !! Whereas MEF90 uses (0,0), (1,0), (0,1), so we need to rescale the affine transform
                  Bt%XX = BBinv(1) * 0.5_kr; Bt%XY = BBinv(3) * 0.5_kr
                  Bt%YX = BBinv(2) * 0.5_kr; Bt%YY = BBinv(4) * 0.5_kr
                  detBBinv = 4.0_kr * abs(detBBinv)
                  call ElementPLagrange2DScalInit(dElem(iELoc), Bt, detBBinv, elemType%order, dQuadratureOrder, ierr)
               end do
               deallocate (BBinv)
               deallocate (BB)
               deallocate (v0)
            end if
            PetscCall(ISRestoreIndices(CellIS, CellID, ierr))
         case (MEF90P1Lagrange2DBoundary%name, MEF90P2Lagrange2DBoundary%name)
            PetscCall(ISGetIndices(CellIS, CellID, ierr))
            allocate (dElem(size(cellID)), stat=ierr)
            if (size(CellID) > 0) then
               allocate (centroid(2))
               allocate (innerNormal(2))
               do iELoc = 1, size(CellID)
                  PetscCall(DMPlexComputeCellGeometryFVM(dm, cellID(iEloc), length, centroid, innerNormal, ierr))
                  outerNormal = -innerNormal
                  call ElementPLagrange2DBoundaryScalInit(dElem(iELoc), length, outerNormal, elemType%order, dQuadratureOrder, ierr)
               end do
               deallocate (innerNormal)
               deallocate (centroid)
            end if
            PetscCall(ISRestoreIndices(CellIS, CellID, ierr))
            !Case (MEF90Q1Lagrange2D%name,MEF90Q2Lagrange2D%name,MEF90Q1Lagrange2DBoundary%name,MEF90Q2Lagrange2DBoundary%name)
            !   !!! Get quadrature points for the current element using DMPlexComputeCellGeometryFEM
            !   !!! Initialize element
         case Default
            write (*, *) __FUNCT__, ': Element type not implemented yet', elemType%name
            ierr = PETSC_ERR_SUP
         end select
      end subroutine Element2DScalInitSet

#undef __FUNCT__
#define __FUNCT__ "Element2DVectInitSet"
      subroutine Element2DVectInitSet(dm, cellIS, dElem, dQuadratureOrder, elemType, ierr)
         type(tDM), intent(IN)                             :: dm
         type(tIS), intent(IN)                             :: cellIS
         type(MEF90Element2DVect), dimension(:), pointer    :: dElem
         PetscInt, intent(IN)                              :: dQuadratureOrder
         type(MEF90ElementType), intent(IN)                :: elemType
         PetscErrorCode, intent(OUT)                       :: ierr

         PetscInt                                         :: iELoc
         PetscInt, dimension(:), pointer                    :: CellID
         PetscReal, dimension(:), pointer                 :: v0, BB, BBinv
         PetscReal                                        :: detBBinv
         type(Mat2D)                                      :: Bt
         PetscReal                                        :: length
         PetscReal, dimension(:), pointer                   :: centroid, innerNormal
         type(Vect2D)                                     :: outerNormal

         select case (elemType%name)
         case (MEF90P1Lagrange2D%name, MEF90P2Lagrange2D%name)
            PetscCall(ISGetIndices(CellIS, CellID, ierr))
            allocate (dElem(size(cellID)), stat=ierr)
            if (size(CellID) > 0) then
               allocate (v0(2))
               allocate (BB(4))
               allocate (BBinv(4))
               Do_Elem_iE: do iELoc = 1, size(CellID)
                  PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cellID(iELoc), v0, BB, BBinv, detBBinv, ierr))
                  !! Petsc uses a reference simplex with vertices at (-1,-1), (1,-1) and (-1,1)
                  !! Whereas MEF90 uses (0,0), (1,0), (0,1), so we need to rescale the affine transform
                  Bt%XX = BBinv(1) * 0.5_kr; Bt%XY = BBinv(3) * 0.5_kr
                  Bt%YX = BBinv(2) * 0.5_kr; Bt%YY = BBinv(4) * 0.5_kr
                  detBBinv = 4.0_kr * abs(detBBinv)
                  call ElementPLagrange2DVectInit(dElem(iELoc), Bt, detBBinv, elemType%order, dQuadratureOrder, ierr)
               end do Do_Elem_iE
               deallocate (BBinv)
               deallocate (BB)
               deallocate (v0)
            end if
            PetscCall(ISRestoreIndices(CellIS, CellID, ierr))
         case (MEF90P1Lagrange2DBoundary%name, MEF90P2Lagrange2DBoundary%name)
            PetscCall(ISGetIndices(CellIS, CellID, ierr))
            allocate (dElem(size(cellID)), stat=ierr)
            if (size(CellID) > 0) then
               allocate (centroid(2))
               allocate (innerNormal(2))
               do iELoc = 1, size(CellID)
                  PetscCall(DMPlexComputeCellGeometryFVM(dm, cellID(iEloc), length, centroid, innerNormal, ierr))
                  outerNormal = -innerNormal
                  call ElementPLagrange2DBoundaryVectInit(dElem(iELoc), length, outerNormal, elemType%order, dQuadratureOrder, ierr)
               end do
               deallocate (innerNormal)
               deallocate (centroid)
            end if
            PetscCall(ISRestoreIndices(CellIS, CellID, ierr))
            !Case (MEF90Q1Lagrange2D%name,MEF90Q2Lagrange2D%name,MEF90Q1Lagrange2DBoundary%name,MEF90Q2Lagrange2DBoundary%name)
            !   !!! Get quadrature points for the current element using DMPlexComputeCellGeometryFEM
            !   !!! Initialize element
         case Default
            write (*, *) __FUNCT__, ': Element type not implemented yet', elemType%name
            ierr = PETSC_ERR_SUP
         end select
      end subroutine Element2DVectInitSet

#undef __FUNCT__
#define __FUNCT__ "Element3DScalInitSet"
      subroutine Element3DScalInitSet(dm, cellIS, dElem, dQuadratureOrder, elemType, ierr)
         type(tDM), intent(IN)                             :: dm
         type(tIS), intent(IN)                             :: cellIS
         type(MEF90Element3DScal), dimension(:), pointer    :: dElem
         PetscInt, intent(IN)                              :: dQuadratureOrder
         type(MEF90ElementType), intent(IN)                :: elemType
         PetscErrorCode, intent(OUT)                       :: ierr

         PetscInt                                         :: iELoc
         PetscInt, dimension(:), pointer                    :: CellID
         PetscReal, dimension(:), pointer                 :: v0, BB, BBinv
         PetscReal                                        :: detBBinv
         type(Mat3D)                                      :: Bt
         PetscReal                                        :: area
         PetscReal, dimension(:), pointer                   :: centroid, innerNormal
         type(Vect3D)                                     :: outerNormal

         select case (elemType%name)
         case (MEF90P1Lagrange3D%name, MEF90P2Lagrange3D%name)
            PetscCall(ISGetIndices(CellIS, CellID, ierr))
            allocate (dElem(size(cellID)), stat=ierr)
            if (size(CellID) > 0) then
               allocate (v0(3))
               allocate (BB(9))
               allocate (BBinv(9))
               Do_Elem_iE: do iELoc = 1, size(CellID)
                  PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cellID(iELoc), v0, BB, BBinv, detBBinv, ierr))
                  !! Petsc uses a reference simplex with vertices at (-1,-1), (1,-1) and (-1,1)
                  !! Whereas MEF90 uses (0,0), (1,0), (0,1), so we need to rescale the affine transform
                  Bt%XX = BBinv(1) * 0.5_kr; Bt%XY = BBinv(4) * 0.5_kr; Bt%XZ = BBinv(7) * 0.5_kr
                  Bt%YX = BBinv(2) * 0.5_kr; Bt%YY = BBinv(5) * 0.5_kr; Bt%YZ = BBinv(8) * 0.5_kr
                  Bt%ZX = BBinv(3) * 0.5_kr; Bt%ZY = BBinv(6) * 0.5_kr; Bt%ZZ = BBinv(9) * 0.5_kr
                  detBBinv = 8.0_kr * abs(detBBinv)
                  call ElementPLagrange3DScalInit(dElem(iELoc), Bt, detBBinv, elemType%order, dQuadratureOrder, ierr)
               end do Do_Elem_iE
               deallocate (BBinv)
               deallocate (BB)
               deallocate (v0)
            end if
            PetscCall(ISRestoreIndices(CellIS, CellID, ierr))
         case (MEF90P1Lagrange3DBoundary%name, MEF90P2Lagrange3DBoundary%name)
            PetscCall(ISGetIndices(CellIS, CellID, ierr))
            allocate (dElem(size(cellID)), stat=ierr)
            if (size(CellID) > 0) then
               allocate (centroid(3))
               allocate (innerNormal(3))
               do iELoc = 1, size(CellID)
                  PetscCall(DMPlexComputeCellGeometryFVM(dm, cellID(iEloc), area, centroid, innerNormal, ierr))
                  outerNormal = -innerNormal
                  call ElementPLagrange3DBoundaryScalInit(dElem(iELoc), area, outerNormal, elemType%order, dQuadratureOrder, ierr)
               end do
               deallocate (innerNormal)
               deallocate (centroid)
            end if
            PetscCall(ISRestoreIndices(CellIS, CellID, ierr))
            !Case (MEF90Q1Lagrange3D%name,MEF90Q2Lagrange3D%name,MEF90Q1Lagrange3DBoundary%name,MEF90Q2Lagrange3DBoundary%name)
            !   !!! Get quadrature points for the current element using DMPlexComputeCellGeometryFEM
            !   !!! Initialize element
         case Default
            write (*, *) __FUNCT__, ': Element type not implemented yet', elemType%name
            ierr = PETSC_ERR_SUP
         end select
      end subroutine Element3DScalInitSet

#undef __FUNCT__
#define __FUNCT__ "Element3DVectInitSet"
      subroutine Element3DVectInitSet(dm, cellIS, dElem, dQuadratureOrder, elemType, ierr)
         type(tDM), intent(IN)                             :: dm
         type(tIS), intent(IN)                             :: cellIS
         type(MEF90Element3DVect), dimension(:), pointer    :: dElem
         PetscInt, intent(IN)                              :: dQuadratureOrder
         type(MEF90ElementType), intent(IN)                :: elemType
         PetscErrorCode, intent(OUT)                       :: ierr

         PetscInt                                         :: iELoc
         PetscInt, dimension(:), pointer                    :: CellID
         PetscReal, dimension(:), pointer                 :: v0, BB, BBinv
         PetscReal                                        :: detBBinv
         type(Mat3D)                                      :: Bt
         PetscReal                                        :: area
         PetscReal, dimension(:), pointer                   :: centroid, innerNormal
         type(Vect3D)                                     :: outerNormal

         select case (elemType%name)
         case (MEF90P1Lagrange3D%name, MEF90P2Lagrange3D%name)
            PetscCall(ISGetIndices(CellIS, CellID, ierr))
            allocate (dElem(size(cellID)), stat=ierr)
            if (size(CellID) > 0) then
               allocate (v0(3))
               allocate (BB(9))
               allocate (BBinv(9))
               Do_Elem_iE: do iELoc = 1, size(CellID)
                  PetscCall(DMPlexComputeCellGeometryAffineFEM(dm, cellID(iELoc), v0, BB, BBinv, detBBinv, ierr))
                  !! Petsc uses a reference simplex with vertices at (-1,-1), (1,-1) and (-1,1)
                  !! Whereas MEF90 uses (0,0), (1,0), (0,1), so we need to rescale the affine transform
                  Bt%XX = BBinv(1) * 0.5_kr; Bt%XY = BBinv(4) * 0.5_kr; Bt%XZ = BBinv(7) * 0.5_kr
                  Bt%YX = BBinv(2) * 0.5_kr; Bt%YY = BBinv(5) * 0.5_kr; Bt%YZ = BBinv(8) * 0.5_kr
                  Bt%ZX = BBinv(3) * 0.5_kr; Bt%ZY = BBinv(6) * 0.5_kr; Bt%ZZ = BBinv(9) * 0.5_kr
                  detBBinv = 8.0_kr * abs(detBBinv)
                  call ElementPLagrange3DVectInit(dElem(iELoc), Bt, detBBinv, elemType%order, dQuadratureOrder, ierr)
               end do Do_Elem_iE
               deallocate (BBinv)
               deallocate (BB)
               deallocate (v0)
            end if
            PetscCall(ISRestoreIndices(CellIS, CellID, ierr))
         case (MEF90P1Lagrange3DBoundary%name, MEF90P2Lagrange3DBoundary%name)
            PetscCall(ISGetIndices(CellIS, CellID, ierr))
            allocate (dElem(size(cellID)), stat=ierr)
            if (size(CellID) > 0) then
               allocate (centroid(3))
               allocate (innerNormal(3))
               do iELoc = 1, size(CellID)
                  PetscCall(DMPlexComputeCellGeometryFVM(dm, cellID(iEloc), area, centroid, innerNormal, ierr))
                  outerNormal = -innerNormal
                  call ElementPLagrange3DBoundaryVectInit(dElem(iELoc), area, outerNormal, elemType%order, dQuadratureOrder, ierr)
               end do
               deallocate (innerNormal)
               deallocate (centroid)
            end if
            PetscCall(ISRestoreIndices(CellIS, CellID, ierr))
            !Case (MEF90Q1Lagrange3D%name,MEF90Q2Lagrange3D%name,MEF90Q1Lagrange3DBoundary%name,MEF90Q2Lagrange3DBoundary%name)
            !   !!! Get quadrature points for the current element using DMPlexComputeCellGeometryFEM
            !   !!! Initialize element
         case Default
            write (*, *) __FUNCT__, ': Element type not implemented yet', elemType%name
            ierr = PETSC_ERR_SUP
         end select
      end subroutine Element3DVectInitSet

#undef __FUNCT__
#define __FUNCT__ "ElementPLagrange2DScalInit"
      subroutine ElementPLagrange2DScalInit(dElem, Bt, DetBinv, dPolynomialOrder, dQuadratureOrder, ierr)
         ! Compute the quadrature weights and the value of the basis functions and their gradient
         ! at the quadrature points.
         ! Quadrature rules courtesy of Shawn Walker walker@math.lsu.edu
         type(MEF90Element2DScal), intent(INOUT) :: dElem
         type(Mat2D), intent(IN)                 :: Bt          ! The transposed of transformation matrix
         PetscReal, intent(IN)                   :: DetBinv     ! The determinant of B^{-1}
         PetscInt                               :: dPolynomialOrder, dQuadratureOrder
         PetscErrorCode, intent(OUT)             :: ierr

         PetscInt                               :: Nb_Gauss
         PetscInt                               :: Num_Dof
         PetscInt                               :: iDoF, iG

         PetscReal, dimension(:, :), pointer       :: PhiHat      ! PhiHat(i,k) The value of the ith basis function at the kth integration point
         type(Vect2D), dimension(:, :), pointer    :: GradPhiHat

         type(Vect2D), dimension(:), pointer     :: Xi ! The quadrature points coordinates in the reference element

         select case (dQuadratureOrder)
         case (0, 1)
            Nb_Gauss = 1
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1)%X = 1.0_kr / 3.0_kr
            Xi(1)%Y = 1.0_kr / 3.0_kr
            dElem%Gauss_C = detBinv / 2.0_kr
         case (2)
            Nb_Gauss = 3
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [1.0_kr / 6.0_kr, 1.0_kr / 6.0_kr]
            Xi(2) = [2.0_kr / 3.0_kr, 1.0_kr / 6.0_kr]
            Xi(3) = [1.0_kr / 6.0_kr, 2.0_kr / 3.0_kr]
            dElem%Gauss_C(1:Nb_Gauss) = detBinv / 6.0_kr
            dElem%Gauss_C(Nb_Gauss) = detBinv / 2.0_kr - sum(dElem%Gauss_C(1:Nb_Gauss - 1))
         case (3)
            Nb_Gauss = 4
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [1.0_kr / 3.0_kr, 1.0_kr / 3.0_kr]
            Xi(2) = [3.0_kr / 5.0_kr, 1.0_kr / 5.0_kr]
            Xi(3) = [1.0_kr / 5.0_kr, 3.0_kr / 5.0_kr]
            Xi(4) = [1.0_kr / 5.0_kr, 1.0_kr / 5.0_kr]
            dElem%Gauss_C(1) = -detBinv * 9.0_kr / 32.0_kr
            dElem%Gauss_C(2:4) = detBinv * 25.0_kr / 96.0_kr
            dElem%Gauss_C(Nb_Gauss) = detBinv / 2.0_kr - sum(dElem%Gauss_C(1:Nb_Gauss - 1))
         case (4)
            Nb_Gauss = 6
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.816847572980459_kr, 0.091576213509771_kr]
            Xi(2) = [0.091576213509771_kr, 0.816847572980459_kr]
            Xi(3) = [0.091576213509771_kr, 0.091576213509771_kr]
            Xi(4) = [0.108103018168070_kr, 0.445948490915965_kr]
            Xi(5) = [0.445948490915965_kr, 0.108103018168070_kr]
            Xi(6) = [0.445948490915965_kr, 0.445948490915965_kr]
            dElem%Gauss_C(1:3) = 0.109951743655322 / 2.0_kr * detBinv
            dElem%Gauss_C(4:6) = 0.223381589678011 / 2.0_kr * detBinv
            dElem%Gauss_C(Nb_Gauss) = detBinv / 2.0_kr - sum(dElem%Gauss_C(1:Nb_Gauss - 1))
         case (5)
            Nb_Gauss = 7
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.3333333333333333_kr, 0.3333333333333333_kr]
            Xi(2) = [0.4701420641051150_kr, 0.4701420641051150_kr]
            Xi(3) = [0.0597158717897698_kr, 0.4701420641051150_kr]
            Xi(4) = [0.4701420641051150_kr, 0.0597158717897698_kr]
            Xi(5) = [0.1012865073234563_kr, 0.1012865073234563_kr]
            Xi(6) = [0.7974269853530872_kr, 0.1012865073234563_kr]
            Xi(7) = [0.1012865073234563_kr, 0.7974269853530872_kr]
            dElem%Gauss_C(1) = .225_kr / 2.0_kr * detBinv
            dElem%Gauss_C(2:4) = 0.13239415278850616 / 2.0_kr * detBinv
            dElem%Gauss_C(5:7) = 0.12593918054482717 / 2.0_kr * detBinv
            dElem%Gauss_C(Nb_Gauss) = detBinv / 2.0_kr - sum(dElem%Gauss_C(1:Nb_Gauss - 1))
         case (-6)
         !! It seems to me that this is really a quadrature rule of order 5...
            Nb_Gauss = 9
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.124949503233232_kr, 0.437525248383384_kr]
            Xi(2) = [0.437525248383384_kr, 0.124949503233232_kr]
            Xi(3) = [0.437525248383384_kr, 0.437525248383384_kr]
            Xi(4) = [0.797112651860071_kr, 0.165409927389841_kr]
            Xi(5) = [0.797112651860071_kr, 0.037477420750088_kr]
            Xi(6) = [0.165409927389841_kr, 0.797112651860071_kr]
            Xi(7) = [0.165409927389841_kr, 0.037477420750088_kr]
            Xi(8) = [0.037477420750088_kr, 0.797112651860071_kr]
            Xi(9) = [0.037477420750088_kr, 0.165409927389841_kr]
            dElem%Gauss_C(1:3) = 0.205950504760887_kr / 2.0_kr * detBinv
            dElem%Gauss_C(4:9) = 0.063691414286223_kr / 2.0_kr * detBinv
            dElem%Gauss_C(Nb_Gauss) = detBinv / 2.0_kr - sum(dElem%Gauss_C(1:Nb_Gauss - 1))
         case (6)
            Nb_Gauss = 12
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.873821971016996_kr, 0.063089014491502_kr]
            Xi(2) = [0.063089014491502_kr, 0.873821971016996_kr]
            Xi(3) = [0.063089014491502_kr, 0.063089014491502_kr]
            Xi(4) = [0.501426509658179_kr, 0.249286745170910_kr]
            Xi(5) = [0.249286745170910_kr, 0.501426509658179_kr]
            Xi(6) = [0.249286745170910_kr, 0.249286745170910_kr]
            Xi(7) = [0.636502499121399_kr, 0.310352451033785_kr]
            Xi(8) = [0.636502499121399_kr, 0.053145049844816_kr]
            Xi(9) = [0.310352451033785_kr, 0.636502499121399_kr]
            Xi(10) = [0.310352451033785_kr, 0.053145049844816_kr]
            Xi(11) = [0.053145049844816_kr, 0.636502499121399_kr]
            Xi(12) = [0.053145049844816_kr, 0.310352451033785_kr]
            dElem%Gauss_C(1:3) = 0.050844906370207_kr / 2.0_kr * detBinv
            dElem%Gauss_C(4:6) = 0.116786275726379_kr / 2.0_kr * detBinv
            dElem%Gauss_C(7:12) = 0.082851075618374_kr / 2.0_kr * detBinv
            dElem%Gauss_C(Nb_Gauss) = detBinv / 2.0_kr - sum(dElem%Gauss_C(1:Nb_Gauss - 1))
         case (7)
            Nb_Gauss = 13
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.333333333333333_kr, 0.333333333333333_kr]
            Xi(2) = [0.479308067841923_kr, 0.260345966079038_kr]
            Xi(3) = [0.260345966079038_kr, 0.479308067841923_kr]
            Xi(4) = [0.260345966079038_kr, 0.260345966079038_kr]
            Xi(5) = [0.869739794195568_kr, 0.065130102902216_kr]
            Xi(6) = [0.065130102902216_kr, 0.869739794195568_kr]
            Xi(7) = [0.065130102902216_kr, 0.065130102902216_kr]
            Xi(8) = [0.638444188569809_kr, 0.312865496004875_kr]
            Xi(9) = [0.638444188569809_kr, 0.048690315425316_kr]
            Xi(10) = [0.312865496004875_kr, 0.638444188569809_kr]
            Xi(11) = [0.312865496004875_kr, 0.048690315425316_kr]
            Xi(12) = [0.048690315425316_kr, 0.638444188569809_kr]
            Xi(13) = [0.048690315425316_kr, 0.312865496004875_kr]
            dElem%Gauss_C(1) = -0.149570044467670_kr / 2.0_kr * detBinv
            dElem%Gauss_C(2:4) = 0.175615257433204_kr / 2.0_kr * detBinv
            dElem%Gauss_C(5:7) = 0.053347235608839_kr / 2.0_kr * detBinv
            dElem%Gauss_C(8:13) = 0.077113760890257_kr / 2.0_kr * detBinv
            dElem%Gauss_C(Nb_Gauss) = detBinv / 2.0_kr - sum(dElem%Gauss_C(1:Nb_Gauss - 1))
         case (8, 9)
            Nb_Gauss = 19
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.3333333333333333_kr, 0.3333333333333333_kr]
            Xi(2) = [0.0206349616025259_kr, 0.4896825191987370_kr]
            Xi(3) = [0.4896825191987370_kr, 0.0206349616025259_kr]
            Xi(4) = [0.4896825191987370_kr, 0.4896825191987370_kr]
            Xi(5) = [0.1258208170141290_kr, 0.4370895914929355_kr]
            Xi(6) = [0.4370895914929355_kr, 0.1258208170141290_kr]
            Xi(7) = [0.4370895914929355_kr, 0.4370895914929355_kr]
            Xi(8) = [0.6235929287619356_kr, 0.1882035356190322_kr]
            Xi(9) = [0.1882035356190322_kr, 0.6235929287619356_kr]
            Xi(10) = [0.1882035356190322_kr, 0.1882035356190322_kr]
            Xi(11) = [0.9105409732110941_kr, 0.0447295133944530_kr]
            Xi(12) = [0.0447295133944530_kr, 0.9105409732110941_kr]
            Xi(13) = [0.0447295133944530_kr, 0.0447295133944530_kr]
            Xi(14) = [0.7411985987844980_kr, 0.0368384120547363_kr]
            Xi(15) = [0.7411985987844980_kr, 0.2219628891607657_kr]
            Xi(16) = [0.0368384120547363_kr, 0.7411985987844980_kr]
            Xi(17) = [0.0368384120547363_kr, 0.2219628891607657_kr]
            Xi(18) = [0.2219628891607657_kr, 0.7411985987844980_kr]
            Xi(19) = [0.2219628891607657_kr, 0.0368384120547363_kr]
            dElem%Gauss_C(1) = 0.09713579628279610_kr / 2.0_kr * detBinv
            dElem%Gauss_C(2:4) = 0.03133470022713983_kr / 2.0_kr * detBinv
            dElem%Gauss_C(5:7) = 0.07782754100477543_kr / 2.0_kr * detBinv
            dElem%Gauss_C(8:10) = 0.07964773892720910_kr / 2.0_kr * detBinv
            dElem%Gauss_C(11:13) = 0.02557767565869810_kr / 2.0_kr * detBinv
            dElem%Gauss_C(14:19) = 0.04328353937728940_kr / 2.0_kr * detBinv
            dElem%Gauss_C(Nb_Gauss) = detBinv / 2.0_kr - sum(dElem%Gauss_C(1:Nb_Gauss - 1))
         case Default
            write (*, *) __FUNCT__, ': Unimplemented quadrature order', dQuadratureOrder
            stop
         end select

         select case (dPolynomialOrder)
         case (1)
            Num_DoF = 3
            allocate (PhiHat(Num_DoF, Nb_Gauss), stat=ierr)
            allocate (GradPhiHat(Num_DoF, Nb_Gauss), stat=ierr)
            PhiHat(1, :) = 1.0_kr - Xi%X - Xi%Y
            PhiHat(2, :) = Xi%X
            PhiHat(3, :) = Xi%Y

            GradPhiHat(1, :)%X = -1.0_kr; GradPhiHat(1, :)%Y = -1.0_kr
            GradPhiHat(2, :)%X = 1.0_kr; GradPhiHat(2, :)%Y = 0.0_kr
            GradPhiHat(3, :)%X = 0.0_kr; GradPhiHat(3, :)%Y = 1.0_kr

         case (2)
            Num_DoF = 6
            allocate (PhiHat(Num_DoF, Nb_Gauss), stat=ierr)
            allocate (GradPhiHat(Num_DoF, Nb_Gauss), stat=ierr)
         !! See Dof local Ordering.md for the unusual dof ordering
            PhiHat(4, :) = (1.0_kr - Xi%X - Xi%Y) * (1.0_kr - 2.0_kr * Xi%X - 2.0_kr * Xi%Y)
            PhiHat(5, :) = Xi%X * (2.0_kr * Xi%X - 1.0_kr)
            PhiHat(6, :) = Xi%Y * (2.0_kr * Xi%Y - 1.0_kr)
            PhiHat(1, :) = 4.0_kr * Xi%X * (1.0_kr - Xi%X - Xi%Y)
            PhiHat(2, :) = 4.0_kr * Xi%X * Xi%Y
            PhiHat(3, :) = 4.0_kr * Xi%Y * (1.0_kr - Xi%X - Xi%Y)

            GradPhiHat(4, :)%X = 4.0_kr * Xi%X + 4.0_kr * Xi%y - 3.0_kr; GradPhiHat(4, :)%Y = 4.0_kr * Xi%X + 4.0_kr * Xi%y - 3.0_kr
            GradPhiHat(5, :)%X = 4.0_kr * Xi%X - 1.0_kr; GradPhiHat(5, :)%Y = 0.0_kr
            GradPhiHat(6, :)%X = 0.0_kr; GradPhiHat(6, :)%Y = 4.0_kr * Xi%Y - 1.0_kr
            GradPhiHat(1, :)%X = 4.0_kr * (1.0_kr - 2.0_kr * Xi%X - Xi%Y); GradPhiHat(1, :)%Y = -4.0_kr * Xi%X
            GradPhiHat(2, :)%X = 4.0_kr * Xi%Y; GradPhiHat(2, :)%Y = 4.0_kr * Xi%X
            GradPhiHat(3, :)%X = -4.0_kr * Xi%Y; GradPhiHat(3, :)%Y = 4.0_kr * (1.0_kr - Xi%X - 2.0_kr * Xi%Y)
         case Default
            Num_DoF = 0
            write (*, *) __FUNCT__, ': Unimplemented PolynomialOrder', dPolynomialOrder
            stop
         end select

         allocate (dElem%BF(Num_DoF, Nb_Gauss), stat=ierr)
         allocate (dElem%Grad_BF(Num_DoF, Nb_Gauss), stat=ierr)
         dElem%BF = PhiHat
         do iDoF = 1, Num_DoF
            do iG = 1, Nb_Gauss
               dElem%Grad_BF(iDoF, iG) = Bt * GradPhiHat(iDoF, iG)
            end do
         end do

         dElem%outerNormal = 0.0_kr

         deallocate (Xi, stat=ierr)
         deallocate (PhiHat, stat=ierr)
         deallocate (GradPhiHat, stat=ierr)
      end subroutine ElementPLagrange2DScalInit

#undef __FUNCT__
#define __FUNCT__ "ElementPLagrange2DBoundaryScalInit"
      subroutine ElementPLagrange2DBoundaryScalInit(dElem, l, outerNormal, dPolynomialOrder, dQuadratureOrder, ierr)
         type(MEF90Element2DScal), intent(INOUT) :: dElem
         PetscReal, intent(IN)                   :: l
         type(Vect2D), intent(IN)                :: outerNormal
         PetscInt, intent(IN)                    :: dPolynomialOrder, dQuadratureOrder
         PetscErrorCode, intent(OUT)             :: ierr

         PetscReal, dimension(:), pointer       :: Xi
         PetscInt                               :: iDoF, iG, Num_Gauss, Num_DoF

         num_Dof = 0
         num_Gauss = 0
         dElem%outerNormal = outerNormal
         select case (dQuadratureOrder)
         case (0, 1)
            Num_Gauss = 1
            allocate (Xi(Num_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Num_Gauss), stat=ierr)
            Xi(1) = 0.0_kr
            dElem%Gauss_C(1) = l
         case (2, 3)
            Num_Gauss = 2
            allocate (Xi(Num_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Num_Gauss), stat=ierr)
            Xi(1) = -1.0_kr / sqrt(3.0_kr)
            Xi(2) = 1.0_kr / sqrt(3.0_kr)
            dElem%Gauss_C = l*.5_kr
         case (4, 5)
            Num_Gauss = 3
            allocate (Xi(Num_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Num_Gauss), stat=ierr)
            Xi(1) = -sqrt(15.0_kr) / 5.0_kr
            Xi(2) = 0.0_kr
            Xi(3) = sqrt(15.0_kr) / 5.0_kr
            dElem%Gauss_C(1) = l * 5.0_kr / 18.0_kr
            dElem%Gauss_C(2) = l * 4.0_kr / 9.0_kr
            dElem%Gauss_C(3) = l * 5.0_kr / 18.0_kr
         case (6, 7)
            Num_Gauss = 4
            allocate (Xi(Num_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Num_Gauss), stat=ierr)
            Xi(1) = -sqrt((15.0_kr + 2.0_kr * sqrt(30.0_kr)) / 35.0_kr)
            Xi(2) = -sqrt((15.0_kr - 2.0_kr * sqrt(30.0_kr)) / 35.0_kr)
            Xi(3) = sqrt((15.0_kr - 2.0_kr * sqrt(30.0_kr)) / 35.0_kr)
            Xi(4) = sqrt((15.0_kr + 2.0_kr * sqrt(30.0_kr)) / 35.0_kr)
            dElem%Gauss_C(1) = l * (.25_kr - sqrt(5.0_kr / 864.0_kr))
            dElem%Gauss_C(2) = l * (.25_kr + sqrt(5.0_kr / 864.0_kr))
            dElem%Gauss_C(3) = l * (.25_kr + sqrt(5.0_kr / 864.0_kr))
            dElem%Gauss_C(4) = l * (.25_kr - sqrt(5.0_kr / 864.0_kr))
         case (8, 9, 10, 11, 12)
            Num_Gauss = 5
            allocate (Xi(Num_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Num_Gauss), stat=ierr)
            Xi(1) = -sqrt(5.0_kr + sqrt(40.0_kr / 7.0_kr)) / 3.0_kr
            Xi(2) = -sqrt(5.0_kr - sqrt(40.0_kr / 7.0_kr)) / 3.0_kr
            Xi(3) = 0.0_kr
            Xi(4) = sqrt(5.0_kr - sqrt(40.0_kr / 7.0_kr)) / 3.0_kr
            Xi(5) = sqrt(5.0_kr + sqrt(40.0_kr / 7.0_kr)) / 3.0_kr
            dElem%Gauss_C(1) = l * (322.0_kr - sqrt(11830.0_kr)) / 1800.0_kr
            dElem%Gauss_C(2) = l * (322.0_kr + sqrt(11830.0_kr)) / 1800.0_kr
            dElem%Gauss_C(3) = l * 128.0_kr / 450.0_kr
            dElem%Gauss_C(4) = l * (322.0_kr + sqrt(11830.0_kr)) / 1800.0_kr
            dElem%Gauss_C(5) = l * (322.0_kr - sqrt(11830.0_kr)) / 1800.0_kr
         case Default
            write (*, *) __FUNCT__, ': Unimplemented quadrature order', dQuadratureOrder
            stop
         end select
         select case (dPolynomialOrder)
         case (1)
            Num_DoF = 2
            allocate (dElem%BF(Num_DoF, Num_Gauss), stat=ierr)
            dElem%BF(1, :) = (1.0_kr - Xi)*.5_kr
            dElem%BF(2, :) = (1.0_kr + Xi)*.5_kr
         case (2)
            Num_DoF = 3
            allocate (dElem%BF(Num_DoF, Num_Gauss), stat=ierr)
         !! See Dof local Ordering.md for the unusual dof ordering
            dElem%BF(2, :) = Xi * (Xi - 1.0_kr)*.5_kr
            dElem%BF(3, :) = Xi * (Xi + 1.0_kr)*.5_kr
            dElem%BF(1, :) = (1.0_kr - Xi) * (1.0_kr + Xi)
         case Default
            write (*, *) '[ERROR]: Polynomial order ', dPolynomialOrder, ' not implemented in ', __FUNCT__
         end select
         allocate (delem%Grad_BF(Num_DoF, Num_Gauss), stat=ierr)
         do iDof = 1, num_dof
            do iG = 1, num_Gauss
               delem%Grad_BF(iDof, iG) = 0.0_kr
            end do
         end do
         deallocate (Xi, stat=ierr)
      end subroutine ElementPLagrange2DBoundaryScalInit

#undef __FUNCT__
#define __FUNCT__ "ElementPLagrange2DVectInit"
      subroutine ElementPLagrange2DVectInit(dElem, Bt, DetBinv, dPolynomialOrder, dQuadratureOrder, ierr)
         type(MEF90Element2DVect), intent(INOUT) :: dElem
         type(Mat2D), intent(IN)                 :: Bt          ! The transposed of transformation matrix
         PetscReal, intent(IN)                   :: DetBinv     ! The determinant of B^{-1}
         PetscInt, intent(IN)                    :: dPolynomialOrder, dQuadratureOrder
         PetscErrorCode, intent(OUT)             :: ierr

         type(MEF90Element2DScal)              :: Elem_Scal
         PetscInt                               :: dim = 2
         PetscInt                               :: Num_DoF, Nb_Gauss, i, iDof, iG

         call ElementPLagrange2DScalInit(Elem_Scal, Bt, DetBinv, dPolynomialOrder, dQuadratureOrder, ierr)
         Num_DoF = size(Elem_Scal%BF, 1)
         Nb_Gauss = size(Elem_Scal%BF, 2)

         allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
         allocate (dElem%BF(Num_DoF * dim, Nb_Gauss), stat=ierr)
         allocate (dElem%Grad_BF(Num_DoF * dim, Nb_Gauss), stat=ierr)
         allocate (dElem%GradS_BF(Num_DoF * dim, Nb_Gauss), stat=ierr)

         dElem%Gauss_C = Elem_Scal%Gauss_C
         do iDof = 1, num_dof * dim
            do iG = 1, Nb_Gauss
               dElem%BF(iDof, iG) = 0.0_kr
               delem%Grad_BF(iDof, iG) = 0.0_kr
               delem%GradS_BF(iDof, iG) = 0.0_kr
            end do
         end do

         do i = 0, Num_DoF - 1
            dElem%BF(i * dim + 1, :)%X = Elem_Scal%BF(i + 1, :)
            dElem%BF(i * dim + 2, :)%Y = Elem_Scal%BF(i + 1, :)

            dElem%Grad_BF(i * dim + 1, :)%XX = Elem_Scal%Grad_BF(i + 1, :)%X
            dElem%Grad_BF(i * dim + 1, :)%XY = Elem_Scal%Grad_BF(i + 1, :)%Y
            dElem%Grad_BF(i * dim + 2, :)%YX = Elem_Scal%Grad_BF(i + 1, :)%X
            dElem%Grad_BF(i * dim + 2, :)%YY = Elem_Scal%Grad_BF(i + 1, :)%Y

            dElem%GradS_BF(i * dim + 1, :)%XX = Elem_Scal%Grad_BF(i + 1, :)%X
            dElem%GradS_BF(i * dim + 1, :)%XY = Elem_Scal%Grad_BF(i + 1, :)%Y / 2.0_kr
            dElem%GradS_BF(i * dim + 2, :)%XY = Elem_Scal%Grad_BF(i + 1, :)%X / 2.0_kr
            dElem%GradS_BF(i * dim + 2, :)%YY = Elem_Scal%Grad_BF(i + 1, :)%Y
         end do

         call MEF90ElementDestroy(Elem_Scal, ierr)
      end subroutine ElementPLagrange2DVectInit

#undef __FUNCT__
#define __FUNCT__ "ElementPLagrange2DBoundaryVectInit"
      subroutine ElementPLagrange2DBoundaryVectInit(dElem, length, outerNormal, dPolynomialOrder, dQuadratureOrder, ierr)
         type(MEF90Element2DVect), intent(INOUT) :: dElem
         PetscReal, intent(IN)                   :: length
         type(Vect2D), intent(IN)                :: outerNormal
         PetscInt, intent(IN)                    :: dPolynomialOrder, dQuadratureOrder
         PetscErrorCode, intent(OUT)             :: ierr

         type(MEF90Element2DScal)              :: Elem_Scal
         PetscInt                               :: dim = 2
         PetscInt                               :: Num_DoF, Nb_Gauss, iDof, iG

         call ElementPLagrange2DBoundaryScalInit(Elem_Scal, length, outerNormal, dPolynomialOrder, dQuadratureOrder, ierr)
         Num_DoF = size(Elem_Scal%BF, 1)
         Nb_Gauss = size(Elem_Scal%BF, 2)
         allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
         allocate (dElem%BF(Num_DoF * dim, Nb_Gauss), stat=ierr)
         allocate (dElem%Grad_BF(Num_DoF * dim, Nb_Gauss), stat=ierr)
         allocate (dElem%GradS_BF(Num_DoF * dim, Nb_Gauss), stat=ierr)

         dElem%Gauss_C = Elem_Scal%Gauss_C
         do iDof = 1, num_dof * dim
            do iG = 1, Nb_Gauss
               dElem%BF(iDof, iG) = 0.0_kr
               delem%Grad_BF(iDof, iG) = 0.0_kr
               delem%GradS_BF(iDof, iG) = 0.0_kr
            end do
         end do

         do idof = 0, Num_DoF - 1
            dElem%BF(iDof * dim + 1, :)%X = Elem_Scal%BF(iDof + 1, :)
            dElem%BF(iDof * dim + 2, :)%Y = Elem_Scal%BF(iDof + 1, :)
         end do

         dElem%outerNormal = Elem_Scal%outerNormal
         call MEF90ElementDestroy(Elem_Scal, ierr)
      end subroutine ElementPLagrange2DBoundaryVectInit

#undef __FUNCT__
#define __FUNCT__ "ElementPLagrange3DScalInit"
      subroutine ElementPLagrange3DScalInit(dElem, Bt, DetBinv, dPolynomialOrder, dQuadratureOrder, ierr)
         ! Compute the quadrature weights and the value of the basis functions and their gradient
         ! at the quadrature points.
         type(MEF90Element3DScal), intent(INOUT) :: dElem
         type(Mat3D), intent(IN)                 :: Bt          ! The transposed of transformation matrix
         PetscReal, intent(IN)                   :: DetBinv     ! The determinant of B^{-1}
         PetscInt, intent(IN)                    :: dPolynomialOrder, dQuadratureOrder
         PetscErrorCode, intent(OUT)             :: ierr

         PetscInt                               :: Nb_Gauss
         PetscInt                               :: Num_Dof
         PetscInt                               :: iDoF, iG

         PetscReal, dimension(:, :), pointer       :: PhiHat      ! PhiHat(i,k) The value of the ith basis function at the kth integration point
         type(Vect3D), dimension(:, :), pointer    :: GradPhiHat

         type(Vect3D), dimension(:), pointer      :: Xi          ! The quadrature points coordinates in the reference element
         PetscReal                              :: a, b         ! Location of integration points in Aiken p. 272 table 10.4

         select case (dQuadratureOrder)
         case (0, 1)
            Nb_Gauss = 1
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [.25_kr, .25_kr, .25_kr]
            dElem%Gauss_C(1) = 1.0_kr / 6.0_kr * detBinv
         case (2)
            a = (5.0_kr + 3.0_kr * sqrt(5.0_kr)) / 20.0_kr
            b = (5.0_kr - sqrt(5.0_kr)) / 20.0_kr
            Nb_Gauss = 4
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [a, b, b]
            Xi(2) = [b, a, b]
            Xi(3) = [b, b, a]
            Xi(4) = [b, b, b]
            dElem%Gauss_C(1:4) = 1.0_kr / 24.0_kr * detBinv
         case (3)
            Nb_Gauss = 5
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [.25_kr, .25_kr, .25_kr]
            Xi(2) = [.5_kr, 1._kr / 6._kr, 1._kr / 6._kr]
            Xi(3) = [1._kr / 6._kr, .5_kr, 1._kr / 6._kr]
            Xi(4) = [1._kr / 6._kr, 1._kr / 6._kr, .5_kr]
            Xi(5) = [1._kr / 6._kr, 1._kr / 6._kr, 1._kr / 6._kr]
            dElem%Gauss_C(1) = -2.0_kr / 15.0_kr * detBinv
            dElem%Gauss_C(2:5) = 3.0_kr / 40.0_kr * detBinv
         case (4)
            Nb_Gauss = 11
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.2500000000000000_kr, 0.2500000000000000_kr, 0.2500000000000000_kr]
            Xi(2) = [0.7857142857142857_kr, 0.0714285714285714_kr, 0.0714285714285714_kr]
            Xi(3) = [0.0714285714285714_kr, 0.0714285714285714_kr, 0.0714285714285714_kr]
            Xi(4) = [0.0714285714285714_kr, 0.0714285714285714_kr, 0.7857142857142857_kr]
            Xi(5) = [0.0714285714285714_kr, 0.7857142857142857_kr, 0.0714285714285714_kr]
            Xi(6) = [0.1005964238332008_kr, 0.3994035761667992_kr, 0.3994035761667992_kr]
            Xi(7) = [0.3994035761667992_kr, 0.1005964238332008_kr, 0.3994035761667992_kr]
            Xi(8) = [0.3994035761667992_kr, 0.3994035761667992_kr, 0.1005964238332008_kr]
            Xi(9) = [0.3994035761667992_kr, 0.1005964238332008_kr, 0.1005964238332008_kr]
            Xi(10) = [0.1005964238332008_kr, 0.3994035761667992_kr, 0.1005964238332008_kr]
            Xi(11) = [0.1005964238332008_kr, 0.1005964238332008_kr, 0.3994035761667992_kr]
            dElem%Gauss_C(1) = -0.0789333333333333_kr / 6.0_kr * detBinv
            dElem%Gauss_C(2:5) = 0.0457333333333333_kr / 6.0_kr * detBinv
            dElem%Gauss_C(6:11) = 0.1493333333333333_kr / 6.0_kr * detBinv
         case (5)
            Nb_Gauss = 15
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.2500000000000000_kr, 0.2500000000000000_kr, 0.2500000000000000_kr]
            Xi(2) = [0.0000000000000000_kr, 0.3333333333333333_kr, 0.3333333333333333_kr]
            Xi(3) = [0.3333333333333333_kr, 0.3333333333333333_kr, 0.3333333333333333_kr]
            Xi(4) = [0.3333333333333333_kr, 0.3333333333333333_kr, 0.0000000000000000_kr]
            Xi(5) = [0.3333333333333333_kr, 0.0000000000000000_kr, 0.3333333333333333_kr]
            Xi(6) = [0.7272727272727273_kr, 0.0909090909090909_kr, 0.0909090909090909_kr]
            Xi(7) = [0.0909090909090909_kr, 0.0909090909090909_kr, 0.0909090909090909_kr]
            Xi(8) = [0.0909090909090909_kr, 0.0909090909090909_kr, 0.7272727272727273_kr]
            Xi(9) = [0.0909090909090909_kr, 0.7272727272727273_kr, 0.0909090909090909_kr]
            Xi(10) = [0.4334498464263357_kr, 0.0665501535736643_kr, 0.0665501535736643_kr]
            Xi(11) = [0.0665501535736643_kr, 0.4334498464263357_kr, 0.0665501535736643_kr]
            Xi(12) = [0.0665501535736643_kr, 0.0665501535736643_kr, 0.4334498464263357_kr]
            Xi(13) = [0.0665501535736643_kr, 0.4334498464263357_kr, 0.4334498464263357_kr]
            Xi(14) = [0.4334498464263357_kr, 0.0665501535736643_kr, 0.4334498464263357_kr]
            Xi(15) = [0.4334498464263357_kr, 0.4334498464263357_kr, 0.0665501535736643_kr]
            dElem%Gauss_C(1) = 0.1817020685825351_kr / 6.0_kr * detBinv
            dElem%Gauss_C(2:5) = 0.0361607142857143_kr / 6.0_kr * detBinv
            dElem%Gauss_C(6:9) = 0.0698714945161738_kr / 6.0_kr * detBinv
            dElem%Gauss_C(10:15) = 0.0656948493683187_kr / 6.0_kr * detBinv
         case (6)
            Nb_Gauss = 24
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.3561913862225449_kr, 0.2146028712591517_kr, 0.2146028712591517_kr]
            Xi(2) = [0.2146028712591517_kr, 0.2146028712591517_kr, 0.2146028712591517_kr]
            Xi(3) = [0.2146028712591517_kr, 0.2146028712591517_kr, 0.3561913862225449_kr]
            Xi(4) = [0.2146028712591517_kr, 0.3561913862225449_kr, 0.2146028712591517_kr]
            Xi(5) = [0.8779781243961660_kr, 0.0406739585346113_kr, 0.0406739585346113_kr]
            Xi(6) = [0.0406739585346113_kr, 0.0406739585346113_kr, 0.0406739585346113_kr]
            Xi(7) = [0.0406739585346113_kr, 0.0406739585346113_kr, 0.8779781243961660_kr]
            Xi(8) = [0.0406739585346113_kr, 0.8779781243961660_kr, 0.0406739585346113_kr]
            Xi(9) = [0.0329863295731731_kr, 0.3223378901422757_kr, 0.3223378901422757_kr]
            Xi(10) = [0.3223378901422757_kr, 0.3223378901422757_kr, 0.3223378901422757_kr]
            Xi(11) = [0.3223378901422757_kr, 0.3223378901422757_kr, 0.0329863295731731_kr]
            Xi(12) = [0.3223378901422757_kr, 0.0329863295731731_kr, 0.3223378901422757_kr]
            Xi(13) = [0.2696723314583159_kr, 0.0636610018750175_kr, 0.0636610018750175_kr]
            Xi(14) = [0.0636610018750175_kr, 0.2696723314583159_kr, 0.0636610018750175_kr]
            Xi(15) = [0.0636610018750175_kr, 0.0636610018750175_kr, 0.2696723314583159_kr]
            Xi(16) = [0.6030056647916491_kr, 0.0636610018750175_kr, 0.0636610018750175_kr]
            Xi(17) = [0.0636610018750175_kr, 0.6030056647916491_kr, 0.0636610018750175_kr]
            Xi(18) = [0.0636610018750175_kr, 0.0636610018750175_kr, 0.6030056647916491_kr]
            Xi(19) = [0.0636610018750175_kr, 0.2696723314583159_kr, 0.6030056647916491_kr]
            Xi(20) = [0.2696723314583159_kr, 0.6030056647916491_kr, 0.0636610018750175_kr]
            Xi(21) = [0.6030056647916491_kr, 0.0636610018750175_kr, 0.2696723314583159_kr]
            Xi(22) = [0.0636610018750175_kr, 0.6030056647916491_kr, 0.2696723314583159_kr]
            Xi(23) = [0.2696723314583159_kr, 0.0636610018750175_kr, 0.6030056647916491_kr]
            Xi(24) = [0.6030056647916491_kr, 0.2696723314583159_kr, 0.0636610018750175_kr]
            dElem%Gauss_C(1:4) = 0.0399227502581679_kr / 6.0_kr * detBinv
            dElem%Gauss_C(5:8) = 0.0100772110553207_kr / 6.0_kr * detBinv
            dElem%Gauss_C(9:12) = 0.0553571815436544_kr / 6.0_kr * detBinv
            dElem%Gauss_C(13:24) = 0.0482142857142857_kr / 6.0_kr * detBinv
         case (7)
            Nb_Gauss = 31
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.2500000000000000_kr, 0.2500000000000000_kr, 0.2500000000000000_kr]
            Xi(2) = [0.7653604230090441_kr, 0.0782131923303186_kr, 0.0782131923303186_kr]
            Xi(3) = [0.0782131923303186_kr, 0.0782131923303186_kr, 0.0782131923303186_kr]
            Xi(4) = [0.0782131923303186_kr, 0.0782131923303186_kr, 0.7653604230090441_kr]
            Xi(5) = [0.0782131923303186_kr, 0.7653604230090441_kr, 0.0782131923303186_kr]
            Xi(6) = [0.6344703500082868_kr, 0.1218432166639044_kr, 0.1218432166639044_kr]
            Xi(7) = [0.1218432166639044_kr, 0.1218432166639044_kr, 0.1218432166639044_kr]
            Xi(8) = [0.1218432166639044_kr, 0.1218432166639044_kr, 0.6344703500082868_kr]
            Xi(9) = [0.1218432166639044_kr, 0.6344703500082868_kr, 0.1218432166639044_kr]
            Xi(10) = [0.0023825066607383_kr, 0.3325391644464206_kr, 0.3325391644464206_kr]
            Xi(11) = [0.3325391644464206_kr, 0.3325391644464206_kr, 0.3325391644464206_kr]
            Xi(12) = [0.3325391644464206_kr, 0.3325391644464206_kr, 0.0023825066607383_kr]
            Xi(13) = [0.3325391644464206_kr, 0.0023825066607383_kr, 0.3325391644464206_kr]
            Xi(14) = [0.0000000000000000_kr, 0.5000000000000000_kr, 0.5000000000000000_kr]
            Xi(15) = [0.5000000000000000_kr, 0.0000000000000000_kr, 0.5000000000000000_kr]
            Xi(16) = [0.5000000000000000_kr, 0.5000000000000000_kr, 0.0000000000000000_kr]
            Xi(17) = [0.5000000000000000_kr, 0.0000000000000000_kr, 0.0000000000000000_kr]
            Xi(18) = [0.0000000000000000_kr, 0.5000000000000000_kr, 0.0000000000000000_kr]
            Xi(19) = [0.0000000000000000_kr, 0.0000000000000000_kr, 0.5000000000000000_kr]
            Xi(20) = [0.2000000000000000_kr, 0.1000000000000000_kr, 0.1000000000000000_kr]
            Xi(21) = [0.1000000000000000_kr, 0.2000000000000000_kr, 0.1000000000000000_kr]
            Xi(22) = [0.1000000000000000_kr, 0.1000000000000000_kr, 0.2000000000000000_kr]
            Xi(23) = [0.6000000000000000_kr, 0.1000000000000000_kr, 0.1000000000000000_kr]
            Xi(24) = [0.1000000000000000_kr, 0.6000000000000000_kr, 0.1000000000000000_kr]
            Xi(25) = [0.1000000000000000_kr, 0.1000000000000000_kr, 0.6000000000000000_kr]
            Xi(26) = [0.1000000000000000_kr, 0.2000000000000000_kr, 0.6000000000000000_kr]
            Xi(27) = [0.2000000000000000_kr, 0.6000000000000000_kr, 0.1000000000000000_kr]
            Xi(28) = [0.6000000000000000_kr, 0.1000000000000000_kr, 0.2000000000000000_kr]
            Xi(29) = [0.1000000000000000_kr, 0.6000000000000000_kr, 0.2000000000000000_kr]
            Xi(30) = [0.2000000000000000_kr, 0.1000000000000000_kr, 0.6000000000000000_kr]
            Xi(31) = [0.6000000000000000_kr, 0.2000000000000000_kr, 0.1000000000000000_kr]
            dElem%Gauss_C(1) = 0.1095853407966528_kr / 6.0_kr * detBinv
            dElem%Gauss_C(2:5) = 0.0635996491464850_kr / 6.0_kr * detBinv
            dElem%Gauss_C(6:9) = -0.3751064406859797_kr / 6.0_kr * detBinv
            dElem%Gauss_C(10:13) = 0.0293485515784412_kr / 6.0_kr * detBinv
            dElem%Gauss_C(14:19) = 0.0058201058201058_kr / 6.0_kr * detBinv
            dElem%Gauss_C(20:31) = 0.1653439153439105_kr / 6.0_kr * detBinv
         case (8)
            Nb_Gauss = 45
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.2500000000000000_kr, 0.2500000000000000_kr, 0.2500000000000000_kr]
            Xi(2) = [0.6175871903000830_kr, 0.1274709365666390_kr, 0.1274709365666390_kr]
            Xi(3) = [0.1274709365666390_kr, 0.1274709365666390_kr, 0.1274709365666390_kr]
            Xi(4) = [0.1274709365666390_kr, 0.1274709365666390_kr, 0.6175871903000830_kr]
            Xi(5) = [0.1274709365666390_kr, 0.6175871903000830_kr, 0.1274709365666390_kr]
            Xi(6) = [0.9037635088221031_kr, 0.0320788303926323_kr, 0.0320788303926323_kr]
            Xi(7) = [0.0320788303926323_kr, 0.0320788303926323_kr, 0.0320788303926323_kr]
            Xi(8) = [0.0320788303926323_kr, 0.0320788303926323_kr, 0.9037635088221031_kr]
            Xi(9) = [0.0320788303926323_kr, 0.9037635088221031_kr, 0.0320788303926323_kr]
            Xi(10) = [0.4502229043567190_kr, 0.0497770956432810_kr, 0.0497770956432810_kr]
            Xi(11) = [0.0497770956432810_kr, 0.4502229043567190_kr, 0.0497770956432810_kr]
            Xi(12) = [0.0497770956432810_kr, 0.0497770956432810_kr, 0.4502229043567190_kr]
            Xi(13) = [0.0497770956432810_kr, 0.4502229043567190_kr, 0.4502229043567190_kr]
            Xi(14) = [0.4502229043567190_kr, 0.0497770956432810_kr, 0.4502229043567190_kr]
            Xi(15) = [0.4502229043567190_kr, 0.4502229043567190_kr, 0.0497770956432810_kr]
            Xi(16) = [0.3162695526014501_kr, 0.1837304473985499_kr, 0.1837304473985499_kr]
            Xi(17) = [0.1837304473985499_kr, 0.3162695526014501_kr, 0.1837304473985499_kr]
            Xi(18) = [0.1837304473985499_kr, 0.1837304473985499_kr, 0.3162695526014501_kr]
            Xi(19) = [0.1837304473985499_kr, 0.3162695526014501_kr, 0.3162695526014501_kr]
            Xi(20) = [0.3162695526014501_kr, 0.1837304473985499_kr, 0.3162695526014501_kr]
            Xi(21) = [0.3162695526014501_kr, 0.3162695526014501_kr, 0.1837304473985499_kr]
            Xi(22) = [0.0229177878448171_kr, 0.2319010893971509_kr, 0.2319010893971509_kr]
            Xi(23) = [0.2319010893971509_kr, 0.0229177878448171_kr, 0.2319010893971509_kr]
            Xi(24) = [0.2319010893971509_kr, 0.2319010893971509_kr, 0.0229177878448171_kr]
            Xi(25) = [0.5132800333608811_kr, 0.2319010893971509_kr, 0.2319010893971509_kr]
            Xi(26) = [0.2319010893971509_kr, 0.5132800333608811_kr, 0.2319010893971509_kr]
            Xi(27) = [0.2319010893971509_kr, 0.2319010893971509_kr, 0.5132800333608811_kr]
            Xi(28) = [0.2319010893971509_kr, 0.0229177878448171_kr, 0.5132800333608811_kr]
            Xi(29) = [0.0229177878448171_kr, 0.5132800333608811_kr, 0.2319010893971509_kr]
            Xi(30) = [0.5132800333608811_kr, 0.2319010893971509_kr, 0.0229177878448171_kr]
            Xi(31) = [0.2319010893971509_kr, 0.5132800333608811_kr, 0.0229177878448171_kr]
            Xi(32) = [0.0229177878448171_kr, 0.2319010893971509_kr, 0.5132800333608811_kr]
            Xi(33) = [0.5132800333608811_kr, 0.0229177878448171_kr, 0.2319010893971509_kr]
            Xi(34) = [0.7303134278075384_kr, 0.0379700484718286_kr, 0.0379700484718286_kr]
            Xi(35) = [0.0379700484718286_kr, 0.7303134278075384_kr, 0.0379700484718286_kr]
            Xi(36) = [0.0379700484718286_kr, 0.0379700484718286_kr, 0.7303134278075384_kr]
            Xi(37) = [0.1937464752488044_kr, 0.0379700484718286_kr, 0.0379700484718286_kr]
            Xi(38) = [0.0379700484718286_kr, 0.1937464752488044_kr, 0.0379700484718286_kr]
            Xi(39) = [0.0379700484718286_kr, 0.0379700484718286_kr, 0.1937464752488044_kr]
            Xi(40) = [0.0379700484718286_kr, 0.7303134278075384_kr, 0.1937464752488044_kr]
            Xi(41) = [0.7303134278075384_kr, 0.1937464752488044_kr, 0.0379700484718286_kr]
            Xi(42) = [0.1937464752488044_kr, 0.0379700484718286_kr, 0.7303134278075384_kr]
            Xi(43) = [0.0379700484718286_kr, 0.1937464752488044_kr, 0.7303134278075384_kr]
            Xi(44) = [0.7303134278075384_kr, 0.0379700484718286_kr, 0.1937464752488044_kr]
            Xi(45) = [0.1937464752488044_kr, 0.7303134278075384_kr, 0.0379700484718286_kr]
            dElem%Gauss_C(1) = -0.2359620398477557_kr / 6.0_kr * detBinv
            dElem%Gauss_C(2:5) = 0.0244878963560562_kr / 6.0_kr * detBinv
            dElem%Gauss_C(6:9) = 0.0039485206398261_kr / 6.0_kr * detBinv
            dElem%Gauss_C(10:15) = 0.0263055529507371_kr / 6.0_kr * detBinv
            dElem%Gauss_C(16:21) = 0.0829803830550589_kr / 6.0_kr * detBinv
            dElem%Gauss_C(22:33) = 0.0254426245481023_kr / 6.0_kr * detBinv
            dElem%Gauss_C(34:45) = 0.0134324384376852_kr / 6.0_kr * detBinv
         case Default
            write (*, *) __FUNCT__, ': Unimplemented quadrature order', dQuadratureOrder
            stop
         end select

         select case (dPolynomialOrder)
         case (1)
            Num_DoF = 4
            allocate (PhiHat(Num_DoF, Nb_Gauss), stat=ierr)
            allocate (GradPhiHat(Num_DoF, Nb_Gauss), stat=ierr)
            PhiHat(1, :) = 1.0_kr - Xi%X - Xi%Y - Xi%Z
            PhiHat(2, :) = Xi(:)%Y
            PhiHat(3, :) = Xi(:)%X
            PhiHat(4, :) = Xi(:)%Z

            GradPhiHat(1, :)%X = -1.0_kr; GradPhiHat(1, :)%Y = -1.0_kr; GradPhiHat(1, :)%Z = -1.0_kr
            GradPhiHat(2, :)%X = 0.0_kr; GradPhiHat(2, :)%Y = 1.0_kr; GradPhiHat(2, :)%Z = 0.0_kr
            GradPhiHat(3, :)%X = 1.0_kr; GradPhiHat(3, :)%Y = 0.0_kr; GradPhiHat(3, :)%Z = 0.0_kr
            GradPhiHat(4, :)%X = 0.0_kr; GradPhiHat(4, :)%Y = 0.0_kr; GradPhiHat(4, :)%Z = 1.0_kr
         case (2)
            Num_Dof = 10
            allocate (PhiHat(Num_DoF, Nb_Gauss), stat=ierr)
            allocate (GradPhiHat(Num_DoF, Nb_Gauss), stat=ierr)
            PhiHat(1, :) = 4.0_kr * (1.0_kr - Xi%X - Xi%Y - Xi%Z) * Xi%Y
            PhiHat(2, :) = 4.0_kr * Xi%X * Xi%Y
            PhiHat(3, :) = 4.0_kr * (1.0_kr - Xi%X - Xi%Y - Xi%Z) * Xi%X
            PhiHat(4, :) = 4.0_kr * (1.0_kr - Xi%X - Xi%Y - Xi%Z) * Xi%Z
            PhiHat(5, :) = 4.0_kr * Xi%Y * Xi%Z
            PhiHat(6, :) = 4.0_kr * Xi%X * Xi%Z
            PhiHat(7, :) = (1.0_kr - Xi%X - Xi%Y - Xi%Z) * (1.0_kr - 2.0_kr * Xi%X - 2.0_kr * Xi%Y - 2.0_kr * Xi%Z)
            PhiHat(8, :) = Xi%Y * (2.0_kr * Xi%Y - 1.0_kr)
            PhiHat(9, :) = Xi%X * (2.0_kr * Xi%X - 1.0_kr)
            PhiHat(10, :) = Xi%Z * (2.0_kr * Xi%Z - 1.0_kr)

            GradPhiHat(1, :)%X = -4.0_kr * Xi%Y
            GradPhiHat(1, :)%Y = 4.0_kr * (1.0_kr - Xi%X - 2.0_kr * Xi%Y - Xi%Z)
            GradPhiHat(1, :)%Z = -4.0_kr * Xi%Y

            GradPhiHat(2, :)%X = 4.0_kr * Xi%Y
            GradPhiHat(2, :)%Y = 4.0_kr * Xi%X
            GradPhiHat(2, :)%Z = 0.0_kr

            GradPhiHat(3, :)%X = 4.0_kr * (1.0_kr - 2.0_kr * Xi%X - Xi%Y - Xi%Z)
            GradPhiHat(3, :)%Y = -4.0_kr * Xi%X
            GradPhiHat(3, :)%Z = -4.0_kr * Xi%X

            GradPhiHat(4, :)%X = -4.0_kr * Xi%Z
            GradPhiHat(4, :)%Y = -4.0_kr * Xi%Z
            GradPhiHat(4, :)%Z = 4.0_kr * (1.0_kr - Xi%X - Xi%Y - 2.0_kr * Xi%Z)

            GradPhiHat(5, :)%X = 0.0_kr
            GradPhiHat(5, :)%Y = 4.0_kr * Xi%Z
            GradPhiHat(5, :)%Z = 4.0_kr * Xi%Y

            GradPhiHat(6, :)%X = 4.0_kr * Xi%Z
            GradPhiHat(6, :)%Y = 0.0_kr
            GradPhiHat(6, :)%Z = 4.0_kr * Xi%X

            GradPhiHat(7, :)%X = 4.0_kr * Xi%X + 4.0_kr * Xi%Y + 4.0_kr * Xi%Z - 3.0_kr
            GradPhiHat(7, :)%Y = 4.0_kr * Xi%X + 4.0_kr * Xi%Y + 4.0_kr * Xi%Z - 3.0_kr
            GradPhiHat(7, :)%Z = 4.0_kr * Xi%X + 4.0_kr * Xi%Y + 4.0_kr * Xi%Z - 3.0_kr

            GradPhiHat(8, :)%X = 0.0_kr
            GradPhiHat(8, :)%Y = 4.0_kr * Xi%Y - 1.0_kr
            GradPhiHat(8, :)%Z = 0.0_kr

            GradPhiHat(9, :)%X = 4.0_kr * Xi%X - 1.0_kr
            GradPhiHat(9, :)%Y = 0.0_kr
            GradPhiHat(9, :)%Z = 0.0_kr

            GradPhiHat(10, :)%X = 0.0_kr
            GradPhiHat(10, :)%Y = 0.0_kr
            GradPhiHat(10, :)%Z = 4.0_kr * Xi%Z - 1.0_kr

         case Default
            Num_DoF = 0
            write (*, *) __FUNCT__, ': Unimplemented PolynomialOrder', dPolynomialOrder
            stop
         end select

         allocate (dElem%BF(Num_DoF, Nb_Gauss), stat=ierr)
         allocate (dElem%Grad_BF(Num_DoF, Nb_Gauss), stat=ierr)
         dElem%BF = PhiHat
         do iDoF = 1, Num_DoF
            do iG = 1, Nb_Gauss
               dElem%Grad_BF(iDoF, iG) = Bt * GradPhiHat(iDoF, iG)
            end do
         end do

         deallocate (Xi, stat=ierr)
         deallocate (PhiHat, stat=ierr)
         deallocate (GradPhiHat, stat=ierr)
      end subroutine ElementPLagrange3DScalInit

#undef __FUNCT__
#define __FUNCT__ "ElementPLagrange3DBoundaryScalInit"
      subroutine ElementPLagrange3DBoundaryScalInit(dElem, area, outerNormal, dPolynomialOrder, dQuadratureOrder, ierr)
         ! Compute the quadrature weights and the value of the basis functions and their gradient
         ! at the quadrature points.
         type(MEF90Element3DScal), intent(INOUT) :: dElem
         PetscReal, intent(IN)                   :: area
         type(Vect3D), intent(IN)                :: outerNormal
         PetscInt, intent(IN)                    :: dPolynomialOrder, dQuadratureOrder
         PetscErrorCode, intent(OUT)             :: ierr

         PetscInt                               :: Nb_Gauss
         PetscInt                               :: Num_Dof
         PetscInt                               :: iDoF, iG

         type(Vect2D), dimension(:), pointer      :: Xi ! The quadrature points coordinates in the reference element

         num_Dof = 0
         nb_Gauss = 0

         dElem%outerNormal = outerNormal

         select case (dQuadratureOrder)
         case (0, 1)
            Nb_Gauss = 1
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1)%X = 1.0_kr / 3.0_kr
            Xi(1)%Y = 1.0_kr / 3.0_kr
            dElem%Gauss_C = area

         case (2)
            Nb_Gauss = 3
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [1.0_kr / 6.0_kr, 1.0_kr / 6.0_kr]
            Xi(2) = [2.0_kr / 3.0_kr, 1.0_kr / 6.0_kr]
            Xi(3) = [1.0_kr / 6.0_kr, 2.0_kr / 3.0_kr]
            dElem%Gauss_C = area / 3.0_kr

         case (3)
            Nb_Gauss = 4
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            dElem%Gauss_C = area * 25.0_kr / 48.0_kr
            dElem%Gauss_C(1) = -area * 9.0_kr / 16.0_kr
            Xi(1) = [1.0_kr / 3.0_kr, 1.0_kr / 3.0_kr]
            Xi(2) = [3.0_kr / 5.0_kr, 1.0_kr / 5.0_kr]
            Xi(3) = [1.0_kr / 5.0_kr, 3.0_kr / 5.0_kr]
            Xi(4) = [1.0_kr / 5.0_kr, 1.0_kr / 5.0_kr]
         case (4)
            Nb_Gauss = 7
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            dElem%Gauss_C(1) = area / 20.0_kr
            dElem%Gauss_C(2) = area * 2.0_kr / 15.0_kr
            dElem%Gauss_C(3) = area / 20.0_kr
            dElem%Gauss_C(4) = area * 2.0_kr / 15.0_kr
            dElem%Gauss_C(5) = area / 20.0_kr
            dElem%Gauss_C(6) = area * 2.0_kr / 15.0_kr
            dElem%Gauss_C(7) = area * 9.0_kr / 20.0_kr
            Xi(1) = [0.0_kr, 0.0_kr]
            Xi(2) = [0.5_kr, 0.0_kr]
            Xi(3) = [1.0_kr, 0.0_kr]
            Xi(4) = [0.5_kr, 0.5_kr]
            Xi(5) = [0.0_kr, 1.0_kr]
            Xi(6) = [0.0_kr, 0.5_kr]
            Xi(7) = [1.0_kr / 3.0_kr, 1.0_kr / 3.0_kr]
         case (5, 6)
            Nb_Gauss = 9
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.124949503233232_kr, 0.437525248383384_kr]
            Xi(2) = [0.437525248383384_kr, 0.124949503233232_kr]
            Xi(3) = [0.437525248383384_kr, 0.437525248383384_kr]
            Xi(4) = [0.797112651860071_kr, 0.165409927389841_kr]
            Xi(5) = [0.797112651860071_kr, 0.037477420750088_kr]
            Xi(6) = [0.165409927389841_kr, 0.797112651860071_kr]
            Xi(7) = [0.165409927389841_kr, 0.037477420750088_kr]
            Xi(8) = [0.037477420750088_kr, 0.797112651860071_kr]
            Xi(9) = [0.037477420750088_kr, 0.165409927389841_kr]
            dElem%Gauss_C(1:3) = 0.205950504760887_kr * area
            dElem%Gauss_C(4:9) = 0.063691414286223_kr * area
         case (7)
            Nb_Gauss = 13
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.333333333333333_kr, 0.333333333333333_kr]
            Xi(2) = [0.479308067841923_kr, 0.260345966079038_kr]
            Xi(3) = [0.260345966079038_kr, 0.479308067841923_kr]
            Xi(4) = [0.260345966079038_kr, 0.260345966079038_kr]
            Xi(5) = [0.869739794195568_kr, 0.065130102902216_kr]
            Xi(6) = [0.065130102902216_kr, 0.869739794195568_kr]
            Xi(7) = [0.065130102902216_kr, 0.065130102902216_kr]
            Xi(8) = [0.638444188569809_kr, 0.312865496004875_kr]
            Xi(9) = [0.638444188569809_kr, 0.048690315425316_kr]
            Xi(10) = [0.312865496004875_kr, 0.638444188569809_kr]
            Xi(11) = [0.312865496004875_kr, 0.048690315425316_kr]
            Xi(12) = [0.048690315425316_kr, 0.638444188569809_kr]
            Xi(13) = [0.048690315425316_kr, 0.312865496004875_kr]
            dElem%Gauss_C(1) = -0.149570044467670_kr * area
            dElem%Gauss_C(2:4) = 0.175615257433204_kr * area
            dElem%Gauss_C(5:7) = 0.053347235608839_kr * area
            dElem%Gauss_C(8:13) = 0.077113760890257_kr * area
         case (8, 9)
            Nb_Gauss = 19
            allocate (Xi(Nb_Gauss), stat=ierr)
            allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
            Xi(1) = [0.3333333333333333_kr, 0.3333333333333333_kr]
            Xi(2) = [0.0206349616025259_kr, 0.4896825191987370_kr]
            Xi(3) = [0.4896825191987370_kr, 0.0206349616025259_kr]
            Xi(4) = [0.4896825191987370_kr, 0.4896825191987370_kr]
            Xi(5) = [0.1258208170141290_kr, 0.4370895914929355_kr]
            Xi(6) = [0.4370895914929355_kr, 0.1258208170141290_kr]
            Xi(7) = [0.4370895914929355_kr, 0.4370895914929355_kr]
            Xi(8) = [0.6235929287619356_kr, 0.1882035356190322_kr]
            Xi(9) = [0.1882035356190322_kr, 0.6235929287619356_kr]
            Xi(10) = [0.1882035356190322_kr, 0.1882035356190322_kr]
            Xi(11) = [0.9105409732110941_kr, 0.0447295133944530_kr]
            Xi(12) = [0.0447295133944530_kr, 0.9105409732110941_kr]
            Xi(13) = [0.0447295133944530_kr, 0.0447295133944530_kr]
            Xi(14) = [0.7411985987844980_kr, 0.0368384120547363_kr]
            Xi(15) = [0.7411985987844980_kr, 0.2219628891607657_kr]
            Xi(16) = [0.0368384120547363_kr, 0.7411985987844980_kr]
            Xi(17) = [0.0368384120547363_kr, 0.2219628891607657_kr]
            Xi(18) = [0.2219628891607657_kr, 0.7411985987844980_kr]
            Xi(19) = [0.2219628891607657_kr, 0.0368384120547363_kr]
            dElem%Gauss_C(1) = 0.09713579628279610_kr * area
            dElem%Gauss_C(2:4) = 0.03133470022713983_kr * area
            dElem%Gauss_C(5:7) = 0.07782754100477543_kr * area
            dElem%Gauss_C(8:10) = 0.07964773892720910_kr * area
            dElem%Gauss_C(11:13) = 0.02557767565869810_kr * area
            dElem%Gauss_C(14:19) = 0.04328353937728940_kr * area
         case Default
            write (*, *) __FUNCT__, ': Unimplemented quadrature order', dQuadratureOrder
            ierr = PETSC_ERR_SUP
            stop
         end select

         select case (dPolynomialOrder)
         case (1)
            Num_DoF = 3
            allocate (dElem%BF(Num_DoF, Nb_Gauss), stat=ierr)
            dElem%BF(1, :) = 1.0_kr - Xi%X - Xi%Y
            dElem%BF(2, :) = Xi(:)%X
            dElem%BF(3, :) = Xi(:)%Y
         case (2)
            Num_DoF = 6
            allocate (dElem%BF(Num_DoF, Nb_Gauss), stat=ierr)
            dElem%BF(5, :) = (1.0_kr - Xi%X - Xi%Y) * (1.0_kr - 2.0_kr * Xi%X - 2.0_kr * Xi%Y)
            dElem%BF(6, :) = Xi%X * (2.0_kr * Xi%X - 1.0_kr)
            dElem%BF(4, :) = Xi%Y * (2.0_kr * Xi%Y - 1.0_kr)
            dElem%BF(2, :) = 4.0_kr * Xi%X * (1.0_kr - Xi%X - Xi%Y)
            dElem%BF(3, :) = 4.0_kr * Xi%X * Xi%Y
            dElem%BF(1, :) = 4.0_kr * Xi%Y * (1.0_kr - Xi%X - Xi%Y)
         case Default
            write (*, *) __FUNCT__, ': Unimplemented PolynomialOrder', dPolynomialOrder
            ierr = PETSC_ERR_SUP
         end select

         allocate (delem%Grad_BF(Num_DoF, Nb_Gauss), stat=ierr)
         do iDof = 1, num_Dof
            do iG = 1, Nb_Gauss
               delem%Grad_BF(iDof, iG) = 0.0_kr
            end do
         end do

         deallocate (Xi, stat=ierr)
      end subroutine ElementPLagrange3DBoundaryScalInit

#undef __FUNCT__
#define __FUNCT__ "ElementPLagrange3DVectInit"
      subroutine ElementPLagrange3DVectInit(dElem, Bt, detBinv, dPolynomialOrder, dQuadratureOrder, ierr)
         type(MEF90Element3DVect), intent(INOUT) :: dElem
         type(Mat3D), intent(IN)                 :: Bt
         PetscReal, intent(IN)                   :: detBinv
         PetscInt, intent(IN)                    :: dPolynomialOrder, dQuadratureOrder
         PetscErrorCode, intent(OUT)             :: ierr

         type(MEF90Element3DScal)              :: Elem_Scal
         PetscInt                               :: dim = 3
         PetscInt                               :: Num_DoF, Nb_Gauss, i, iDof, iG

         call ElementPLagrange3DScalInit(Elem_Scal, Bt, detBinv, dPolynomialOrder, dQuadratureOrder, ierr)
         Num_DoF = size(Elem_Scal%BF, 1)
         Nb_Gauss = size(Elem_Scal%BF, 2)
         allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
         allocate (dElem%BF(Num_DoF * dim, Nb_Gauss), stat=ierr)
         allocate (dElem%Grad_BF(Num_DoF * dim, Nb_Gauss), stat=ierr)
         allocate (dElem%GradS_BF(Num_DoF * dim, Nb_Gauss), stat=ierr)

         dElem%Gauss_C = Elem_Scal%Gauss_C
         do iDof = 1, num_dof * dim
            do iG = 1, Nb_Gauss
               dElem%BF(iDof, iG) = 0.0_kr
               delem%Grad_BF(iDof, iG) = 0.0_kr
               delem%GradS_BF(iDof, iG) = 0.0_kr
            end do
         end do

         do i = 0, Num_DoF - 1
            dElem%BF(i * dim + 1, :)%X = Elem_Scal%BF(i + 1, :)
            dElem%BF(i * dim + 2, :)%Y = Elem_Scal%BF(i + 1, :)
            dElem%BF(i * dim + 3, :)%Z = Elem_Scal%BF(i + 1, :)
            dElem%Grad_BF(i * dim + 1, :)%XX = Elem_Scal%Grad_BF(i + 1, :)%X
            dElem%Grad_BF(i * dim + 1, :)%XY = Elem_Scal%Grad_BF(i + 1, :)%Y
            dElem%Grad_BF(i * dim + 1, :)%XZ = Elem_Scal%Grad_BF(i + 1, :)%Z
            dElem%Grad_BF(i * dim + 2, :)%YX = Elem_Scal%Grad_BF(i + 1, :)%X
            dElem%Grad_BF(i * dim + 2, :)%YY = Elem_Scal%Grad_BF(i + 1, :)%Y
            dElem%Grad_BF(i * dim + 2, :)%YZ = Elem_Scal%Grad_BF(i + 1, :)%Z
            dElem%Grad_BF(i * dim + 3, :)%ZX = Elem_Scal%Grad_BF(i + 1, :)%X
            dElem%Grad_BF(i * dim + 3, :)%ZY = Elem_Scal%Grad_BF(i + 1, :)%Y
            dElem%Grad_BF(i * dim + 3, :)%ZZ = Elem_Scal%Grad_BF(i + 1, :)%Z

            dElem%GradS_BF(i * dim + 1, :)%XX = Elem_Scal%Grad_BF(i + 1, :)%X
            dElem%GradS_BF(i * dim + 1, :)%XY = Elem_Scal%Grad_BF(i + 1, :)%Y * 0.5_kr
            dElem%GradS_BF(i * dim + 1, :)%XZ = Elem_Scal%Grad_BF(i + 1, :)%Z * 0.5_kr

            dElem%GradS_BF(i * dim + 2, :)%XY = Elem_Scal%Grad_BF(i + 1, :)%X * 0.5_kr
            dElem%GradS_BF(i * dim + 2, :)%YY = Elem_Scal%Grad_BF(i + 1, :)%Y
            dElem%GradS_BF(i * dim + 2, :)%YZ = Elem_Scal%Grad_BF(i + 1, :)%Z * 0.5_kr

            dElem%GradS_BF(i * dim + 3, :)%XZ = Elem_Scal%Grad_BF(i + 1, :)%X * 0.5_kr
            dElem%GradS_BF(i * dim + 3, :)%YZ = Elem_Scal%Grad_BF(i + 1, :)%Y * 0.5_kr
            dElem%GradS_BF(i * dim + 3, :)%ZZ = Elem_Scal%Grad_BF(i + 1, :)%Z
         end do
         call MEF90ElementDestroy(Elem_Scal, ierr)
      end subroutine ElementPLagrange3DVectInit

#undef __FUNCT__
#define __FUNCT__ "ElementPLagrange3DBoundaryVectInit"
      subroutine ElementPLagrange3DBoundaryVectInit(dElem, area, outerNormal, dPolynomialOrder, dQuadratureOrder, ierr)
         type(MEF90Element3DVect), intent(INOUT) :: dElem
         PetscReal, intent(IN)                   :: area
         type(Vect3D), intent(IN)                :: outerNormal
         PetscInt                               :: dPolynomialOrder, dQuadratureOrder
         PetscErrorCode, intent(OUT)             :: ierr

         type(MEF90Element3DScal)              :: Elem_Scal
         PetscInt                               :: dim = 3
         PetscInt                               :: Num_DoF, Nb_Gauss, iDof, iG

         call ElementPLagrange3DBoundaryScalInit(Elem_Scal, area, outerNormal, dPolynomialOrder, dQuadratureOrder, ierr)
         Num_DoF = size(Elem_Scal%BF, 1)
         Nb_Gauss = size(Elem_Scal%BF, 2)
         allocate (dElem%Gauss_C(Nb_Gauss), stat=ierr)
         allocate (dElem%BF(Num_DoF * dim, Nb_Gauss), stat=ierr)
         allocate (dElem%Grad_BF(Num_DoF * dim, Nb_Gauss), stat=ierr)
         allocate (dElem%GradS_BF(Num_DoF * dim, Nb_Gauss), stat=ierr)

         dElem%Gauss_C = Elem_Scal%Gauss_C
         do iDof = 1, num_dof * dim
            do iG = 1, Nb_Gauss
               dElem%BF(iDof, iG) = 0.0_kr
               delem%Grad_BF(iDof, iG) = 0.0_kr
               delem%GradS_BF(iDof, iG) = 0.0_kr
            end do
         end do
         do iDof = 0, Num_DoF - 1
            dElem%BF(iDof * dim + 1, :)%X = Elem_Scal%BF(iDof + 1, :)
            dElem%BF(iDof * dim + 2, :)%Y = Elem_Scal%BF(iDof + 1, :)
            dElem%BF(iDof * dim + 3, :)%Z = Elem_Scal%BF(iDof + 1, :)
         end do

         dElem%outerNormal = Elem_Scal%outerNormal
         call MEF90ElementDestroy(Elem_Scal, ierr)
      end subroutine ElementPLagrange3DBoundaryVectInit

#undef __FUNCT__
#define __FUNCT__ "Element2DScalDestroy"
      subroutine Element2DScalDestroy(dElem, ierr)
         type(MEF90Element2DScal), intent(INOUT) :: dElem
         PetscErrorCode, intent(OUT)             :: ierr

         if (associated(dElem%BF)) then
            deallocate (dElem%BF, stat=ierr)
         end if
         if (associated(dElem%Grad_BF)) then
            deallocate (dElem%Grad_BF, stat=ierr)
         end if
         if (associated(dElem%Gauss_C)) then
            deallocate (dElem%Gauss_C, stat=ierr)
         end if
      end subroutine Element2DScalDestroy

#undef __FUNCT__
#define __FUNCT__ "Element2DVectDestroy"
      subroutine Element2DVectDestroy(dElem, ierr)
         type(MEF90Element2DVect), intent(INOUT) :: dElem
         PetscErrorCode, intent(OUT)             :: ierr

         if (associated(dElem%BF)) then
            deallocate (dElem%BF, stat=ierr)
         end if
         if (associated(dElem%Grad_BF)) then
            deallocate (dElem%Grad_BF, stat=ierr)
         end if
         if (associated(dElem%GradS_BF)) then
            deallocate (dElem%GradS_BF, stat=ierr)
         end if
         if (associated(dElem%Gauss_C)) then
            deallocate (dElem%Gauss_C, stat=ierr)
         end if
      end subroutine Element2DVectDestroy

#undef __FUNCT__
#define __FUNCT__ "Element3DScalDestroy"
      subroutine Element3DScalDestroy(dElem, ierr)
         type(MEF90Element3DScal), intent(INOUT) :: dElem
         PetscErrorCode, intent(OUT)             :: ierr

         if (associated(dElem%BF)) then
            deallocate (dElem%BF, stat=ierr)
         end if
         if (associated(dElem%Grad_BF)) then
            deallocate (dElem%Grad_BF, stat=ierr)
         end if
         if (associated(dElem%Gauss_C)) then
            deallocate (dElem%Gauss_C, stat=ierr)
         end if
      end subroutine Element3DScalDestroy

#undef __FUNCT__
#define __FUNCT__ "Element3DVectDestroy"
      subroutine Element3DVectDestroy(dElem, ierr)
         type(MEF90Element3DVect), intent(INOUT) :: dElem
         PetscErrorCode, intent(OUT)             :: ierr
         if (associated(dElem%BF)) then
            deallocate (dElem%BF, stat=ierr)
         end if
         if (associated(dElem%Grad_BF)) then
            deallocate (dElem%Grad_BF, stat=ierr)
         end if
         if (associated(dElem%GradS_BF)) then
            deallocate (dElem%GradS_BF, stat=ierr)
         end if
         if (associated(dElem%Gauss_C)) then
            deallocate (dElem%Gauss_C, stat=ierr)
         end if
      end subroutine Element3DVectDestroy

#undef __FUNCT__
#define __FUNCT__ "Element2DScalDestroySet"
      subroutine Element2DScalDestroySet(dElem, ierr)
         type(MEF90Element2DScal), dimension(:), pointer    :: dElem
         PetscErrorCode, intent(OUT)                        :: ierr

         PetscInt                                          :: cell

         do cell = 1, size(dElem)
            call MEF90ElementDestroy(dElem(cell), ierr)
         end do
         deallocate (dElem, stat=ierr)
      end subroutine Element2DScalDestroySet

#undef __FUNCT__
#define __FUNCT__ "Element2DVectDestroySet"
      subroutine Element2DVectDestroySet(dElem, ierr)
         type(MEF90Element2DVect), dimension(:), pointer    :: dElem
         PetscErrorCode, intent(OUT)                        :: ierr

         PetscInt                                          :: cell

         do cell = 1, size(dElem)
            call MEF90ElementDestroy(dElem(cell), ierr)
         end do
         deallocate (dElem, stat=ierr)
      end subroutine Element2DVectDestroySet

#undef __FUNCT__
#define __FUNCT__ "Element3DScalDestroySet"
      subroutine Element3DScalDestroySet(dElem, ierr)
         type(MEF90Element3DScal), dimension(:), pointer    :: dElem
         PetscErrorCode, intent(OUT)                        :: ierr

         PetscInt                                          :: cell

         do cell = 1, size(dElem)
            call MEF90ElementDestroy(dElem(cell), ierr)
         end do
         deallocate (dElem, stat=ierr)
      end subroutine Element3DScalDestroySet

#undef __FUNCT__
#define __FUNCT__ "Element3DVectDestroySet"
      subroutine Element3DVectDestroySet(dElem, ierr)
         type(MEF90Element3DVect), dimension(:), pointer    :: dElem
         PetscErrorCode, intent(OUT)                        :: ierr

         PetscInt                                          :: cell

         do cell = 1, size(dElem)
            call MEF90ElementDestroy(dElem(cell), ierr)
         end do
         deallocate (dElem, stat=ierr)
      end subroutine Element3DVectDestroySet

   end module m_MEF90_Elements
