!/*****************************************************************************/
! *
! *  Elmer, A Finite Element Software for Multiphysical Problems
! *
! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
! * 
! * This library is free software; you can redistribute it and/or
! * modify it under the terms of the GNU Lesser General Public
! * License as published by the Free Software Foundation; either
! * version 2.1 of the License, or (at your option) any later version.
! *
! * This library is distributed in the hope that it will be useful,
! * but WITHOUT ANY WARRANTY; without even the implied warranty of
! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
! * Lesser General Public License for more details.
! * 
! * You should have received a copy of the GNU Lesser General Public
! * License along with this library (in file ../LGPL-2.1); if not, write 
! * to the Free Software Foundation, Inc., 51 Franklin Street, 
! * Fifth Floor, Boston, MA  02110-1301  USA
! *
! *****************************************************************************/
!
!/******************************************************************************
! *
! *  Authors: Mikko Byckling, Juha Ruokolainen
! *  Email:   Juha.Ruokolainen@csc.fi
! *  Web:     http://www.csc.fi/elmer
! *  Address: CSC - IT Center for Science Ltd.
! *           Keilaranta 14
! *           02101 Espoo, Finland 
! *
! *  Original Date: 23 Aug 2004
! *
! *****************************************************************************/

!> \ingroup ElmerLib
!> \{

!-------------------------------------------------------------------------------
!>  Module defining mappings for p elements. These include nodal points 
!>  contained by faces and edges, element boundary maps (edges for 2d elements,
!>  faces for 3d) and mappings from faces to edge numbers. Mappings defined in 
!>  this module are compatible with basis functions defined in module 
!>  PElementBase.
!-------------------------------------------------------------------------------

MODULE PElementMaps
  USE Types
  Use GeneralUtils, ONLY : I2S

  IMPLICIT NONE

  ! Private mappings. For access use get[Element][Type]Map(i)
  PRIVATE QuadEdgeMap, TriangleEdgeMap, &
       TetraEdgeMap1, TetraFaceMap1, TetraFaceEdgeMap1, &
       TetraEdgeMap2, TetraFaceMap2, TetraFaceEdgeMap2, &
       BrickEdgeMap, BrickFaceMap, BrickFaceEdgeMap, &
       WedgeEdgeMap, WedgeFaceMap, WedgeFaceEdgeMap, &
       PyramidEdgeMap, PyramidFaceMap, PyramidFaceEdgeMap, &
       MInit, TriangleFaceMap, QuadFaceMap
  ! Mappings
  INTEGER, TARGET, SAVE :: QuadEdgeMap(4,2), TriangleEdgeMap(3,2), &
       TetraEdgeMap1(6,2), TetraFaceMap1(4,3), TetraFaceEdgeMap1(4,3), &
       TetraEdgeMap2(6,2), TetraFaceMap2(4,3), TetraFaceEdgeMap2(4,3),&
       BrickEdgeMap(12,2), BrickFaceMap(6,4), BrickFaceEdgeMap(6,4), &
       WedgeEdgeMap(9,2), WedgeFaceMap(5,4), WedgeFaceEdgeMap(5,4), &
       PyramidEdgeMap(8,2), PyramidFaceMap(5,4), PyramidFaceEdgeMap(5,4), &
       TriangleFaceMap(1,3), QuadFaceMap(1,4), LineEdgeMap(1,2)

  LOGICAL, SAVE :: MInit = .FALSE.
  !$OMP THREADPRIVATE(MInit, QuadEdgeMap, TriangleEdgeMap, LineEdgeMap, &
  !$OMP&              TetraEdgeMap1, TetraFaceMap1, TetraFaceEdgeMap1, &
  !$OMP&              TetraEdgeMap2, TetraFaceMap2, TetraFaceEdgeMap2,&
  !$OMP&              BrickEdgeMap, BrickFaceMap, BrickFaceEdgeMap, &
  !$OMP&              WedgeEdgeMap, WedgeFaceMap, WedgeFaceEdgeMap, &
  !$OMP&              PyramidEdgeMap, PyramidFaceMap, PyramidFaceEdgeMap, TriangleFaceMap, QuadFaceMap)
CONTAINS

  ! MAPPINGS

  ! First some direct mappings to elements. These should not be used directly 
  ! unless element type is implicitly known from context. Better way is to use
  ! getElement[Boundary,Edge,Face]Map -routines.

    ! Call: localEdge = getQuadEdge(i)
    !
    ! Function returns mapping from edge number to edge endpoints 

    FUNCTION getLineEdgeMap(i) RESULT(localEdge)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(2) :: localEdge

      IF (.NOT. MInit) CALL InitializeMappings()
      
      localEdge(:) = LineEdgeMap(i,:)
    END FUNCTION getLineEdgeMap

    ! Call: localEdge = getQuadEdge(i)
    !
    ! Function returns mapping from edge number to edge endpoints 

    FUNCTION getQuadEdgeMap(i) RESULT(localEdge)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(2) :: localEdge

      IF (.NOT. MInit) CALL InitializeMappings()
      
      localEdge(:) = QuadEdgeMap(i,:)
    END FUNCTION getQuadEdgeMap

    ! Call: localFace = getQuadFaceMap(i)
    ! 
    ! Function returns mapping from face number to face nodes
    FUNCTION getQuadFaceMap(i) RESULT(localFace)
      IMPLICIT NONE
      
      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(4) :: localFace

      IF (.NOT. MInit) CALL InitializeMappings()

      localFace(:) = QuadFaceMap(1,:)
    END FUNCTION getQuadFaceMap

    ! Call: localEdge = getTriangleEdge(i)
    ! 
    ! Function returns mapping from edge number to edge endpoints

    FUNCTION getTriangleEdgeMap(i) RESULT(localEdge)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(2) :: localEdge

      IF (.NOT. MInit) CALL InitializeMappings()

      localEdge(:)=TriangleEdgeMap(i,:)
    END FUNCTION getTriangleEdgeMap

    ! Call: localFace = geTriangleFaceMap(i)
    ! 
    ! Function returns mapping from face number to face nodes
    FUNCTION getTriangleFaceMap(i) RESULT(localFace)
      IMPLICIT NONE
      
      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(3) :: localFace

      IF (.NOT. MInit) CALL InitializeMappings()

      localFace(:) = TriangleFaceMap(1,:)
    END FUNCTION getTriangleFaceMap
    
    ! Call: localEdge = getBrickEdgeMap(i)
    ! 
    ! Function returns mapping from edge number to edge endpoints

    FUNCTION getBrickEdgeMap(i) RESULT(localEdge)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(2) :: localEdge

      IF (.NOT. MInit) CALL InitializeMappings()

      localEdge(:) = BrickEdgeMap(i,:)
    END FUNCTION getBrickEdgeMap
    
    ! Call: localFace = getBrickFaceMap(i)
    ! 
    ! Function returns mapping from face number to face nodes

    FUNCTION getBrickFaceMap(i) RESULT(localFace)
      IMPLICIT NONE
      
      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(4) :: localFace

      IF (.NOT. MInit) CALL InitializeMappings()

      localFace(:) = BrickFaceMap(i,:)
    END FUNCTION getBrickFaceMap

    ! Call: localEdge = getFaceEdgeMap(face, localNode)
    !
    ! getFaceEdgeMap returns number of local edge when given face and 
    ! its local node number. Node number is treated as edges beginning point

    FUNCTION getBrickFaceEdgeMap(face, localNode) RESULT(localEdge)
      IMPLICIT NONE
      CHARACTER(:), ALLOCATABLE :: msg

      ! Parameters 
      INTEGER, INTENT(IN) :: face, localNode
      ! Variables
      INTEGER :: localEdge

      IF (.NOT. MInit) CALL InitializeMappings()

      localEdge = BrickFaceEdgeMap(face,localNode)

      IF (localEdge == 0) THEN
         msg = 'Unknown combination node for (face,node)'//I2S(face)//I2S(localNode)
         CALL Fatal('PElementMaps::getBrickFaceEdgeMap', msg)
      END IF
    END FUNCTION getBrickFaceEdgeMap


    FUNCTION getTetraEdgeMap(i,TYPE) RESULT(edge)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: i
      INTEGER, INTENT(IN), OPTIONAL :: TYPE
      INTEGER :: t
      INTEGER, DIMENSION(2) :: edge

      IF (.NOT. MInit) CALL InitializeMappings()

      ! If type not present use default (1)
      t = 1
      IF (PRESENT(TYPE)) t = TYPE

      ! Select edge map by tetra type
      SELECT CASE (t)
      CASE (1)
         edge(:) = TetraEdgeMap1(i,:)
      CASE (2)
         edge(:) = TetraEdgeMap2(i,:)
      CASE DEFAULT
         CALL Fatal('PElementMaps::getTetraEdgeMap','Unknown tetra type')
      END SELECT
    END FUNCTION getTetraEdgeMap


    FUNCTION getTetraFaceMap(i,TYPE) RESULT(face)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: i
      INTEGER, INTENT(IN), OPTIONAL :: TYPE
      INTEGER :: t
      INTEGER, DIMENSION(3) :: face

      IF (.NOT. MInit) CALL InitializeMappings()
      
      ! If type not present use default (1)
      t = 1
      IF (PRESENT(TYPE)) t = TYPE

      ! Select face map by tetra type 
      SELECT CASE(t)
      CASE (1)
         face(:) = TetraFaceMap1(i,:)
      CASE (2)
         face(:) = TetraFaceMap2(i,:)
      CASE DEFAULT 
         CALL Fatal('PElementMaps::getTetraFaceMap','Unknown tetra type')
      END SELECT
    END FUNCTION getTetraFaceMap

    FUNCTION getWedgeEdgeMap(i) RESULT(edge)
      IMPLICIT NONE
      
      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(2) :: edge

      IF (.NOT. MInit) CALL InitializeMappings()

      edge(:) = WedgeEdgeMap(i,:)
    END FUNCTION getWedgeEdgeMap


    FUNCTION getWedgeFaceMap(i) RESULT(face)
      IMPLICIT NONE
      
      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(4) :: face

      IF (.NOT. MInit) CALL InitializeMappings()

      face(:) = WedgeFaceMap(i,:)
    END FUNCTION getWedgeFaceMap


    FUNCTION getPyramidEdgeMap(i) RESULT(edge)
      IMPLICIT NONE

      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(2) :: edge

      IF (.NOT. MInit) CALL InitializeMappings()

      edge(:) = PyramidEdgeMap(i,:)
    END FUNCTION getPyramidEdgeMap


    FUNCTION getPyramidFaceMap(i) RESULT(face)
      IMPLICIT NONE
      
      INTEGER, INTENT(IN) :: i
      INTEGER, DIMENSION(4) :: face

      IF (.NOT. MInit) CALL InitializeMappings()

      face(:) = PyramidFaceMap(i,:)
    END FUNCTION getPyramidFaceMap


!------------------------------------------------------------------------------
!>     Mapping from element local edge or face number to nodes contained in 
!>     that edge or face. 
!------------------------------------------------------------------------------
    FUNCTION getElementBoundaryMap(Element, i) RESULT(map)
!------------------------------------------------------------------------------
!
!  ARGUMENTS:
!    Type(Element_t) :: Element
!      INPUT: Element to get map for
!
!    INTEGER, INTENT(IN) :: i
!      INPUT: Local number of elements edge or face
!
!  FUNCTION VALUE:
!    INTEGER :: map(4)
!       Map containing local node numbers of given local edge or face
!    
!------------------------------------------------------------------------------
      IMPLICIT NONE

      TYPE(Element_t) :: Element
      INTEGER, INTENT(IN) :: i
      
      INTEGER :: map(4)

      IF (.NOT. MInit) CALL InitializeMappings()
      map = 0

      ! Function is not defined for non p elements
      !IF (.NOT. ASSOCIATED(Element % PDefs)) THEN
      !   CALL Warn('PElementMaps::getElementBoundaryMap','Element not p element') 
      !   RETURN
      !END IF

      SELECT CASE(Element % TYPE % ElementCode / 100)
      CASE (3)
         map(1:2) = getTriangleEdgeMap(i)
      CASE (4)
         map(1:2) = getQuadEdgeMap(i)
      CASE (5)
         map(1:3) = getTetraFaceMap(i,Element % PDefs % TetraType)
      CASE (6)
         map(1:4) = getPyramidFaceMap(i)
      CASE (7)
         map(1:4) = getWedgeFaceMap(i)
      CASE (8)
         map(1:4) = getBrickFaceMap(i)
      CASE DEFAULT
         CALL Fatal('PElementMaps::getElementBoundaryMap','Unsupported element type')
      END SELECT
    END FUNCTION getElementBoundaryMap


!------------------------------------------------------------------------------
!>     Mapping from element local face to local edges contained in face. Given
!>     element and local face number this routine returns numbers of local edges
!>     on face. 
!------------------------------------------------------------------------------
    FUNCTION getFaceEdgeMap( Element, i) RESULT(map)
!------------------------------------------------------------------------------
!
!  ARGUMENTS:
!    Type(Element_t) :: Element
!      INPUT: Element to get map for
!
!    INTEGER, INTENT(IN) :: i
!      INPUT: Local number of element face
!
!  FUNCTION VALUE:
!    INTEGER :: map(4)
!       Map containing local numbers of edges on face
!    
!------------------------------------------------------------------------------
      IMPLICIT NONE

      TYPE(Element_t) :: Element
      INTEGER, INTENT(IN) :: i

      INTEGER :: elementCode, map(4)

      elementCode = Element % TYPE % ElementCode

      IF (.NOT. MInit) CALL InitializeMappings()

      ! Function is not defined for non p elements
      !IF (.NOT. ASSOCIATED(Element % PDefs)) THEN
      !   CALL Warn('PElementMaps::getFaceEdgeMap','Element not p element') 
      !   map = 0
      !   RETURN
      !END IF

      SELECT CASE(elementCode / 100)
      CASE (5)
         map = 0
         SELECT CASE (Element % PDefs % TetraType)
         CASE (1)
            map(1:3) = TetraFaceEdgeMap1(i,:)
         CASE (2)
            map(1:3) = TetraFaceEdgeMap2(i,:)
         CASE DEFAULT
            CALL Fatal('PElementMaps::getFaceEdgeMap','Unknown tetra type')
         END SELECT
      CASE (6)
         map(1:4) = PyramidFaceEdgeMap(i,:)
      CASE (7)
         map(1:4) = WedgeFaceEdgeMap(i,:)
      CASE (8)
         map(1:4) = BrickFaceEdgeMap(i,:)
      CASE DEFAULT
         CALL Fatal('PElementMaps::getFaceEdgeMap','Unsupported element type')
      END SELECT
    END FUNCTION getFaceEdgeMap


!------------------------------------------------------------------------------
!>     Get mappings for given element to element edges and their nodes. Given 
!>     element, this routine returns a map containing nodes (endpoints) of
!>     elements edges. 
!------------------------------------------------------------------------------
    SUBROUTINE GetElementEdgeMap( Element, map )
!------------------------------------------------------------------------------
!
!  ARGUMENTS:
!    Type(Element_t) :: Element
!      INPUT: Element to get map for
!
!    INTEGER :: map(:,:)
!       OUTPUT: Map containing local node numbers of local edges
!    
!------------------------------------------------------------------------------
      IMPLICIT NONE
      TYPE(Element_t) :: Element
      INTEGER,  POINTER :: map(:,:)

      IF (.NOT. MInit) CALL InitializeMappings()

      ! Function is not defined for non p elements
      IF (.NOT. ASSOCIATED(Element % PDefs)) THEN
         CALL Warn('PElementMaps::GetElementEdgeMap','Element not p element') 
         map = 0
         RETURN
      END IF

      SELECT CASE (Element % TYPE % ElementCode / 100)
      CASE (2)
         map => LineEdgeMap
      CASE (3)
         map => TriangleEdgeMap
      CASE (4)
         map => QuadEdgeMap
      CASE (5)
         SELECT CASE( Element % PDefs % TetraType )
         CASE (1)
            map => TetraEdgeMap1
         CASE (2)
            map => TetraEdgeMap2
         CASE DEFAULT
            CALL Fatal('PElementMaps::GetElementEdgeMap','Unknown tetra type for p element')
         END SELECT
      CASE (6)
         map => PyramidEdgeMap
      CASE (7)
         map => WedgeEdgeMap
      CASE (8)
         map => BrickEdgeMap
      CASE DEFAULT
         CALL Fatal('PElementMaps::GetElementEdgeMap','Unsupported element type')
      END SELECT
    END SUBROUTINE GetElementEdgeMap
   

!------------------------------------------------------------------------------
!>     Get mappings for given element to element faces and their nodes. Given 
!>     element, this routine returns a map containing nodes (endpoints) of
!>     elements face. 
!------------------------------------------------------------------------------
    SUBROUTINE GetElementFaceMap( Element, faceMap )
!------------------------------------------------------------------------------
!
!  ARGUMENTS:
!    Type(Element_t) :: Element
!      INPUT: Element to get map for
!
!    INTEGER :: map(:,:)
!       OUTPUT: Map containing local node numbers of local faces
!    
!------------------------------------------------------------------------------
      IMPLICIT NONE
      
      TYPE(Element_t) :: Element
      INTEGER, POINTER :: faceMap(:,:)

      IF (.NOT. MInit) CALL InitializeMappings()

      ! Function is not defined for non p elements
      IF (.NOT. ASSOCIATED(Element % PDefs)) THEN
         CALL Warn('PElementMaps::GetElementFaceMap','Element not p element') 
         NULLIFY(faceMap)
         RETURN
      END IF

      SELECT CASE (Element % TYPE % ElementCode / 100)
       CASE (3)
          facemap => TriangleFaceMap
       CASE (4)
          facemap => QuadFaceMap
      CASE (5)
         SELECT CASE( Element % PDefs % TetraType )
         CASE (1)
            faceMap => TetraFaceMap1
         CASE (2)
            faceMap => TetraFaceMap2
         CASE DEFAULT
            CALL Fatal('PElementMaps::GetElementFaceMap','Unknown tetra type for p element')
         END SELECT
      CASE (6)
         faceMap => PyramidFaceMap
      CASE (7)
         faceMap => WedgeFaceMap
      CASE (8)
         faceMap => BrickFaceMap
      CASE DEFAULT
         CALL Fatal('PElementMaps::GetElementFaceMap','Unsupported element type')
      END SELECT
    END SUBROUTINE GetElementFaceMap


!------------------------------------------------------------------------------    
!>     Get mappings for given element to elements faces and their edge. Given 
!>     element, this routine returns a map containing local edge numbers of
!>     elements faces. 
!------------------------------------------------------------------------------    
    SUBROUTINE GetElementFaceEdgeMap( Element, faceEdgeMap )
!------------------------------------------------------------------------------    
!
!  ARGUMENTS:
!    Type(Element_t) :: Element
!      INPUT: Element to get map for
!
!    INTEGER :: map(:,:)
!       OUTPUT: Map containing local edge numbers of local faces
!    
!------------------------------------------------------------------------------
      IMPLICIT NONE

      TYPE(Element_t) :: Element
      INTEGER, POINTER :: faceEdgeMap(:,:)
      
      IF (.NOT. MInit) CALL InitializeMappings()

      ! Function is not defined for non p elements
      IF (.NOT. ASSOCIATED(Element % PDefs)) THEN
         CALL Warn('PElementMaps::GetElementFaceEdgeMap','Element not p element') 
         NULLIFY(faceEdgeMap)
         RETURN
      END IF

      SELECT CASE (Element % TYPE % ElementCode / 100)
      CASE (5)
         SELECT CASE( Element % PDefs % TetraType )
         CASE (1)
            faceEdgeMap => TetraFaceEdgeMap1
         CASE (2)
            faceEdgeMap => TetraFaceEdgeMap2
         CASE DEFAULT
            CALL Fatal('PElementMaps::GetElementFaceEdgeMap','Unknown tetra type for p element')
         END SELECT
      CASE (6)
         faceEdgeMap => PyramidFaceEdgeMap
      CASE (7)
         faceEdgeMap => WedgeFaceEdgeMap
      CASE (8)
         faceEdgeMap => BrickFaceEdgeMap
      CASE DEFAULT
         CALL Fatal('PElementMaps::GetElementFaceEdgeMap','Unsupported element type')
      END SELECT
    END SUBROUTINE GetElementFaceEdgeMap


!------------------------------------------------------------------------------
!>   This subroutine initializes element mappings.
!------------------------------------------------------------------------------
    SUBROUTINE InitializeMappings() 
!------------------------------------------------------------------------------
      IMPLICIT NONE
      
      CALL Info('PElementMaps::InitializeMappings','Initializing mappings for elements',Level=10)

      LineEdgeMap(1,:) = [1,2]

      ! Quad edge mappings
      QuadEdgeMap(1,:) = (/ 1,2 /)
      QuadEdgeMap(2,:) = (/ 2,3 /)
      QuadEdgeMap(3,:) = (/ 4,3 /)
      QuadEdgeMap(4,:) = (/ 1,4 /)

      QuadFaceMap(1,:) = (/ 1,2,3,4 /)

      ! Triangle edge mappings
      TriangleEdgeMap(1,:) = (/ 1,2 /)
      TriangleEdgeMap(2,:) = (/ 2,3 /)
      TriangleEdgeMap(3,:) = (/ 3,1 /)

      TriangleFaceMap(1,:) = (/ 1,2,3 /)

      ! Brick edge mappings
      BrickEdgeMap(1,:) = (/ 1,2 /)
      BrickEdgeMap(2,:) = (/ 2,3 /)
      BrickEdgeMap(3,:) = (/ 4,3 /)
      BrickEdgeMap(4,:) = (/ 1,4 /)
      BrickEdgeMap(5,:) = (/ 5,6 /)
      BrickEdgeMap(6,:) = (/ 6,7 /)
      BrickEdgeMap(7,:) = (/ 8,7 /)
      BrickEdgeMap(8,:) = (/ 5,8 /)
      BrickEdgeMap(9,:) = (/ 1,5 /)
      BrickEdgeMap(10,:) = (/ 2,6 /)
      BrickEdgeMap(11,:) = (/ 3,7 /)
      BrickEdgeMap(12,:) = (/ 4,8 /)

      ! Brick face mappings
      BrickFaceMap(1,:) = (/ 1,2,3,4 /) ! xi,eta
      BrickFaceMap(2,:) = (/ 5,6,7,8 /) ! xi,eta
      BrickFaceMap(3,:) = (/ 1,2,6,5 /) ! xi,zeta
      BrickFaceMap(4,:) = (/ 2,3,7,6 /) ! eta,zeta
      ! BrickFaceMap(5,:) = (/ 3,4,8,7 /)
      BrickFaceMap(5,:) = (/ 4,3,7,8 /)
      ! BrickFaceMap(6,:) = (/ 4,1,5,8 /) 
      BrickFaceMap(6,:) = (/ 1,4,8,5 /)

      BrickFaceEdgeMap(1,:) = (/ 1,2,3,4 /)
      BrickFaceEdgeMap(2,:) = (/ 5,6,7,8 /)    
      BrickFaceEdgeMap(3,:) = (/ 1,10,5,9 /)
      BrickFaceEdgeMap(4,:) = (/ 2,11,6,10 /)
      ! BrickFaceEdgeMap(5,:) = (/ 3,12,7,11 /)
      BrickFaceEdgeMap(5,:) = (/ 3,11,7,12 /)
      ! BrickFaceEdgeMap(6,:) = (/ 4,9,8,12 /)
      BrickFaceEdgeMap(6,:) = (/ 4,12,8,9 /)

      ! Tetra edge mappings (not needed for enforcing parity!)
      ! Type 1
      TetraEdgeMap1(1,:) = (/ 1,2 /)
      TetraEdgeMap1(2,:) = (/ 2,3 /)
      TetraEdgeMap1(3,:) = (/ 1,3 /)
      TetraEdgeMap1(4,:) = (/ 1,4 /)
      TetraEdgeMap1(5,:) = (/ 2,4 /)
      TetraEdgeMap1(6,:) = (/ 3,4 /)
      ! Type 2
      TetraEdgeMap2(1,:) = (/ 1,2 /)
      TetraEdgeMap2(2,:) = (/ 3,2 /)
      TetraEdgeMap2(3,:) = (/ 1,3 /)
      TetraEdgeMap2(4,:) = (/ 1,4 /)
      TetraEdgeMap2(5,:) = (/ 2,4 /)
      TetraEdgeMap2(6,:) = (/ 3,4 /)

      ! Tetra face mappings (not needed for enforcing parity!)
      ! Type 1
      TetraFaceMap1(1,:) = (/ 1,2,3 /)
      TetraFaceMap1(2,:) = (/ 1,2,4 /)
      TetraFaceMap1(3,:) = (/ 2,3,4 /)
      TetraFaceMap1(4,:) = (/ 1,3,4 /)
      ! Type 2 
      TetraFaceMap2(1,:) = (/ 1,3,2 /)
      TetraFaceMap2(2,:) = (/ 1,2,4 /)
      TetraFaceMap2(3,:) = (/ 3,2,4 /)
      TetraFaceMap2(4,:) = (/ 1,3,4 /)

      ! Type 1 
      TetraFaceEdgeMap1(1,:) = (/ 1,2,3 /)
      TetraFaceEdgeMap1(2,:) = (/ 1,5,4 /)
      TetraFaceEdgeMap1(3,:) = (/ 2,6,5 /)
      TetraFaceEdgeMap1(4,:) = (/ 3,6,4 /)
      ! Type 2 
      TetraFaceEdgeMap2(1,:) = (/ 3,2,1 /)
      TetraFaceEdgeMap2(2,:) = (/ 1,5,4 /)
      TetraFaceEdgeMap2(3,:) = (/ 2,5,6 /)
      TetraFaceEdgeMap2(4,:) = (/ 3,6,4 /)

      ! Wedge edge mappings
      WedgeEdgeMap(1,:) = (/ 1,2 /)
      WedgeEdgeMap(2,:) = (/ 2,3 /)
      WedgeEdgeMap(3,:) = (/ 3,1 /)
      WedgeEdgeMap(4,:) = (/ 4,5 /)
      WedgeEdgeMap(5,:) = (/ 5,6 /)
      WedgeEdgeMap(6,:) = (/ 6,4 /)
      WedgeEdgeMap(7,:) = (/ 1,4 /)
      WedgeEdgeMap(8,:) = (/ 2,5 /)
      WedgeEdgeMap(9,:) = (/ 3,6 /)

      ! Wedge face mappings
      WedgeFaceMap(1,:) = (/ 1,2,3,0 /)
      WedgeFaceMap(2,:) = (/ 4,5,6,0 /)
      WedgeFaceMap(3,:) = (/ 1,2,5,4 /)
      WedgeFaceMap(4,:) = (/ 2,3,6,5 /)
      WedgeFaceMap(5,:) = (/ 3,1,4,6 /)

      WedgeFaceEdgeMap(1,:) = (/ 1,2,3,0 /)
      WedgeFaceEdgeMap(2,:) = (/ 4,5,6,0 /)
      WedgeFaceEdgeMap(3,:) = (/ 1,8,4,7 /)
      WedgeFaceEdgeMap(4,:) = (/ 2,9,5,8 /)
      WedgeFaceEdgeMap(5,:) = (/ 3,7,6,9 /)
      
      ! Pyramid edge mappings 
      PyramidEdgeMap(1,:) = (/ 1,2 /)
      PyramidEdgeMap(2,:) = (/ 2,3 /)
      PyramidEdgeMap(3,:) = (/ 4,3 /)
      PyramidEdgeMap(4,:) = (/ 1,4 /)
      PyramidEdgeMap(5,:) = (/ 1,5 /)
      PyramidEdgeMap(6,:) = (/ 2,5 /)
      PyramidEdgeMap(7,:) = (/ 3,5 /)
      PyramidEdgeMap(8,:) = (/ 4,5 /)

      ! Pyramid face mappings
      PyramidFaceMap(1,:) = (/ 1,2,3,4 /)
      PyramidFaceMap(2,:) = (/ 1,2,5,0 /)
      PyramidFaceMap(3,:) = (/ 2,3,5,0 /)
      PyramidFaceMap(4,:) = (/ 3,4,5,0 /)
      PyramidFaceMap(5,:) = (/ 4,1,5,0 /)

      PyramidFaceEdgeMap(1,:) = (/ 1,2,3,4 /)
      PyramidFaceEdgeMap(2,:) = (/ 1,6,5,0 /)
      PyramidFaceEdgeMap(3,:) = (/ 2,7,6,0 /)
      PyramidFaceEdgeMap(4,:) = (/ 3,8,7,0 /)
      PyramidFaceEdgeMap(5,:) = (/ 4,5,8,0 /)

      MInit = .TRUE.
    END SUBROUTINE InitializeMappings

!------------------------------------------------------------------------------
  FUNCTION getEdgeDOFs( Element, p ) RESULT(EdgeDOFs)
!------------------------------------------------------------------------------
    IMPLICIT NONE
!------------------------------------------------------------------------------
    TYPE(Element_t) :: Element
    INTEGER :: EdgeDOFs
    INTEGER, INTENT(IN) :: p

    IF (.NOT. ASSOCIATED(Element % PDefs) ) THEN
       EdgeDOFs = 0
       RETURN
    END IF

    EdgeDOFs = MAX(0, p-1)
!------------------------------------------------------------------------------
  END FUNCTION getEdgeDOFs
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!>     Based on element face polynomial degree p, return degrees of freedom for
!>     given face. 
!------------------------------------------------------------------------------
  FUNCTION getFaceDOFs(Element, p, faceNumber, Face ) RESULT(faceDOFs)
!------------------------------------------------------------------------------
!
!  ARGUMENTS:
!    Type(Element_t), POINTER :: Element
!      INPUT: Element to get face dofs to 
!
!    INTEGER :: p
!      INPUT: Face polynomial degree p
!
!    INTEGER :: faceNumber
!      INPUT: Local number of face for element (important for wedges and 
!        pyramids).
!
!  FUNCTION VALUE:
!    REAL(KIND=dp) :: faceDOFs
!       number of face dofs for Element
!    
!------------------------------------------------------------------------------
    IMPLICIT NONE
    
    TYPE(Element_t) :: Element
    TYPE(Element_t), OPTIONAL :: Face
    INTEGER, INTENT(IN) :: p
    INTEGER, INTENT(IN), OPTIONAL :: faceNumber
    INTEGER :: faceDOFs
    LOGICAL :: SerendipityPBasis

    ! This function is not defined for non p elements
    IF (.NOT. ASSOCIATED(Element % PDefs) ) THEN
       faceDOFs = 0
       RETURN
    END IF
    SerendipityPBasis = Element % PDefs % Serendipity

    faceDOFs = 0
    IF(p<=1) RETURN

    SELECT CASE(Element % TYPE % ElementCode / 100)
    ! Triangle
    CASE (3)
       faceDOFs = (p-1)*(p-2)/2
    ! Quad
    CASE (4)
       IF ( SerendipityPBasis ) THEN
         faceDOFs = (p-2)*(p-3)/2
       ELSE
         faceDOFs = (p-1)**2 ! (p-1)*p/2
       END IF
    ! Tetrahedron
    CASE (5)
       faceDOFs = (p-1)*(p-2)/2
    ! Pyramid
    CASE (6)
       SELECT CASE(faceNumber)
        CASE (1)
          faceDOFs = (p-1)**2 ! (p-1)*p/2
        CASE (2:5)
          faceDOFs = (p-1)*(p-2)/2
       END SELECT
    ! Wedge
    CASE (7)
       SELECT CASE(faceNumber)
       CASE (1,2)
          faceDOFs = (p-1)*(p-2)/2
       CASE (3:5)
          IF(SerendipityPBasis) THEN
            faceDOFs = (p-2)*(p-3)/2
          ELSE
            faceDOFs = (p-1)**2 ! (p-1)*p/2
          END IF
       END SELECT
    ! Brick   
    CASE (8)
       IF(SerendipityPBasis) THEN
         faceDOFs = (p-2)*(p-3)/2
       ELSE
         faceDOFs = (p-1)**2 ! (p-1)*p/2
       END IF
    CASE DEFAULT
      WRITE(Message,'(A,I0)') 'Unsupported p element type: ',Element % TYPE % ElementCode
      CALL Warn('PElementMaps::getFaceDOFs',Message)
      faceDOFs = p
    END SELECT

    faceDOFs = MAX(0, faceDOFs)
  END FUNCTION getFaceDOFs


!------------------------------------------------------------------------------
!> Based on the polynomial degree p of the element, return the number of
!> bubble functions (the count of bubble DOFs). 
!> NOTE: The returned value is not the bubble count for an approximation
!> based on the space Q_p of polynomials of degree at most p in each variable 
!> separately.
!------------------------------------------------------------------------------
  FUNCTION getBubbleDOFs( Element, p) RESULT(bubbleDOFs)
!------------------------------------------------------------------------------
!
!  ARGUMENTS:
!    Type(Element_t), POINTER :: Element
!      INPUT: Element to get bubble dofs to 
!
!    INTEGER :: p
!      INPUT: Element polynomial degree p
!
!  FUNCTION VALUE:
!    REAL(KIND=dp) :: bubbleDOFs
!       number of bubble dofs for Element
!    
!------------------------------------------------------------------------------
    IMPLICIT NONE
    
    TYPE(Element_t) :: Element
    INTEGER, INTENT(IN) :: p
    INTEGER :: bubbleDOFs, i
    LOGICAL :: SerendipityPBasis
    
    ! This function is not defined for non p elements
    IF (.NOT. ASSOCIATED(Element % PDefs) ) THEN
       bubbleDOFs = 0
       RETURN
    END IF

    SerendipityPBasis = Element % PDefs % Serendipity

    ! Select by element type
    bubbleDOFs = 0
    IF(p<=1) RETURN

    SELECT CASE (Element % TYPE % ElementCode / 100)
    ! Line 
    CASE (2)
      BubbleDOFs = p-1
    ! Triangle
    CASE (3)
      BubbleDOFs = (p-2)*(p-1)/2
    ! Quad
    CASE (4)
      IF(SerendipityPBasis) THEN
        BubbleDOFs = (p-2)*(p-3)/2
      ELSE
        BubbleDOFs = (p-1)**2 ! (p-1)*p/2
      END IF
    ! Tetrahedron
    CASE (5)
      BubbleDOFs = (p-3)*(p-2)*(p-1)/6
    ! Pyramid
    CASE (6)
       BubbleDOFs = (p-2)*(p-1)*p/6
    ! Prism
    CASE (7)
      IF(SerendipityPBasis) THEN
        bubbleDOFs = (p-2)*(p-3)*(p-4)/6
      ELSE
        BubbleDOFs = (p-1)**2*(p-2)/2
      END IF
    ! Hexa
    CASE (8)
      IF(SerendipityPBasis) THEN
        bubbleDOFs = (p-3)*(p-4)*(p-5)/6
      ELSE
        BubbleDOFs = (p-1)**3
      END IF
    CASE DEFAULT
       CALL Warn('PElementMaps::getBubbleDOFs','Unsupported p element type')
       bubbleDOFs = p
    END SELECT

    bubbleDOFs = MAX(0, bubbleDOFs)
!------------------------------------------------------------------------------
  END FUNCTION getBubbleDOFs
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!> Checks if given element is a p-element active in a particular solver.   
!------------------------------------------------------------------------------
  FUNCTION isActivePElement(Element,USolver) RESULT(retVal)
!------------------------------------------------------------------------------
    IMPLICIT NONE

    TYPE(Element_t), INTENT(IN) :: Element
    TYPE(Solver_t), POINTER, OPTIONAL :: USolver
    LOGICAL :: retVal

    INTEGER :: m
    TYPE(Solver_t), POINTER :: pSolver
        
    retVal = isPelement(Element)

    ! The solver can have active p-element only if we have p-elements!
    IF(.NOT. retVal) RETURN
    
    IF( PRESENT( USolver ) ) THEN
      pSolver => USolver
    ELSE
      pSolver => CurrentModel % Solver
    END IF
    
    IF(ASSOCIATED(pSolver))THEN
      IF(ALLOCATED(pSolver % Def_Dofs)) THEN
        m = Element % Type % ElementCode / 100
        retVal = ANY(pSolver % Def_Dofs(m,:,6)>0)
      END IF
    END IF

!------------------------------------------------------------------------------
  END FUNCTION isActivePElement
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!> Checks whether given solver has active p-element definitions (for some
!> element type and for some body)
!------------------------------------------------------------------------------
  FUNCTION isActivePSolver(Solver) RESULT(retVal)
!------------------------------------------------------------------------------
    IMPLICIT NONE

    TYPE(Solver_t) :: Solver
    LOGICAL :: retVal
 
    TYPE(Element_t) :: Element
        
    retVal = .FALSE.    
    IF(ALLOCATED(Solver % Def_Dofs)) THEN
      retVal = ANY(Solver % Def_Dofs(:,:,6)>0)
    END IF
    
!------------------------------------------------------------------------------
  END FUNCTION isActivePSolver
!------------------------------------------------------------------------------

  

!------------------------------------------------------------------------------
!> Checks whether given element has p-element information associated  
!------------------------------------------------------------------------------
    FUNCTION isPElement( Element ) RESULT(retVal)
!------------------------------------------------------------------------------
!
!  ARGUMENTS:
!    Type(Element_t) :: Element
!      INPUT: Element to check
!
!  FUNCTION VALUE:
!    LOGICAL :: retVal
!       .TRUE. if given element is a p element, .FALSE. otherwise
!    
!------------------------------------------------------------------------------
      IMPLICIT NONE

      TYPE(Element_t), INTENT(IN) :: Element
      LOGICAL :: retVal

      retVal = ASSOCIATED(Element % PDefs)
!------------------------------------------------------------------------------
    END FUNCTION isPElement
!------------------------------------------------------------------------------
    

!------------------------------------------------------------------------------
  FUNCTION getEffectiveBubbleP(Element,set_p,bdofs) RESULT(p)
!------------------------------------------------------------------------------
     IMPLICIT NONE

     INTEGER :: p, bdofs, set_p
     TYPE(Element_t) :: Element

     INTEGER :: curr_nb, nb

     p = set_p
     curr_nb = GetBubbleDOFs(Element,p)
     nb = MAX( curr_nb, bdofs )

     DO WHILE(curr_nb < nb)
       p = p + 1
       curr_nb = GetBubbleDOFs(Element,p)
     END DO
!------------------------------------------------------------------------------
  END FUNCTION getEffectiveBubbleP
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!> Get the number of Gauss points for P-elements.
!------------------------------------------------------------------------------
    FUNCTION getNumberOfGaussPoints( Element, Mesh ) RESULT(ngp)
!------------------------------------------------------------------------------
      IMPLICIT NONE
      TYPE(Mesh_t) :: Mesh 
      TYPE(Element_t) :: Element
      INTEGER :: ngp
!------------------------------------------------------------------------------
      INTEGER :: edgeP, faceP, bubbleP, TrueBubbleP, maxp

      IF (.NOT. ASSOCIATED(Element % PDefs)) THEN
         CALL Warn('PElementMaps::getNumberOfGaussPoints','Element not p element')
         ngp = 0
         RETURN
      END IF

      ! Max p of edges
      edgeP = 0
      IF ( Element % TYPE % DIMENSION == 2 .OR. &
           Element % TYPE % DIMENSION == 3) THEN
         edgeP = getEdgeP( Element, Mesh )
      END IF
      
      ! Max p of faces
      faceP = 0
      IF ( Element % TYPE % DIMENSION == 3 ) THEN
         faceP = getFaceP(Element, Mesh )
      END IF
      
      ! Element bubble p
      bubbleP = getEffectiveBubbleP(Element,Element % PDefs % P,Element % bdofs)
      TrueBubbleP = bubbleP

      ! Special quadrature may be available: 
      IF (Element % PDefs % Serendipity .AND. Element % TYPE % ElementCode / 100 == 4) THEN
        ! The true polynomial degree is as follows
        !maxp = MAX(1, edgeP, faceP, TrueBubbleP)
        ! but this would replace the true bubble degree by a tampered value:
        maxp = MAX(1, edgeP, faceP, BubbleP)

        ! Economic quadratures cannot be used if an explicit bubble augmentation is used
        ! with lower-order finite elements:
        IF ( .NOT.(Element % PDefs % P<4 .AND. Element % BDOFs>0) ) THEN
          IF (maxp > 1 .AND. maxp <= 7) THEN
            !PRINT *, 'SETTING SPECIAL NGP'
            !PRINT *, 'MAXP=',MAXP
            SELECT CASE(maxp)
            CASE(2)
              ngp = 8
            CASE(3)
              ngp = 12
            CASE(4)
              ngp = 20
            CASE(5)
              ngp = 36
            CASE(6)
              ngp = 45
            CASE(7)
              ngp = 60
            END SELECT
            RETURN
          END IF
        END IF
      END IF
      ! Get the number r of Gauss points for the product of two basis functions: 
      ! r = (2*max(p)+1)/2
      maxp = MAX(1, edgeP, faceP, bubbleP) + 1
      ! The number of Gauss points based on the Cartesian product (more efficient
      ! rules could be defined):
      ngp = maxp ** Element % TYPE % DIMENSION
!------------------------------------------------------------------------------
    END FUNCTION getNumberOfGaussPoints
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
    FUNCTION getEdgeP( Element, Mesh ) RESULT(edgeP)
!------------------------------------------------------------------------------
      IMPLICIT NONE
      
      TYPE(Mesh_t) :: Mesh
      TYPE(Element_t) :: Element 
      TYPE(Element_t), POINTER :: Edge
      
      INTEGER :: edgeP, i

      IF (.NOT. ASSOCIATED(Element % PDefs)) THEN
         CALL Warn('PElementMaps::getEdgeP','Element not p element')
         edgeP = 0
         RETURN
      END IF

      ! Get max p of edges of element if any
      edgeP = 0
      IF (ASSOCIATED(Element % EdgeIndexes)) THEN
         DO i=1, Element % TYPE % NumberOfEdges
            Edge => Mesh % Edges(Element % EdgeIndexes(i))
            ! Here if edge has no dofs it effectively has degree of 0
            IF (Edge % BDOFs <= 0) CYCLE
            edgeP = MAX(edgeP,Edge % PDefs % P)
         END DO
      END IF
!------------------------------------------------------------------------------
    END FUNCTION getEdgeP
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
    FUNCTION getFaceP( Element, Mesh ) RESULT(faceP)
!------------------------------------------------------------------------------
      IMPLICIT NONE

      TYPE(Element_t) :: Element
      TYPE(Element_t), POINTER :: Face
      INTEGER :: faceP, i
      TYPE(Mesh_t) :: Mesh
      
      IF (.NOT. ASSOCIATED(Element % PDefs)) THEN
         CALL Warn('PElementMaps::getFaceP','Element not p element')
         faceP = 0
         RETURN
      END IF

      faceP = 0
      IF (ASSOCIATED(Element % FaceIndexes)) THEN
         DO i=1,Element % TYPE % NumberOfFaces
            Face => Mesh % Faces( Element % FaceIndexes(i) )
            ! Here if face has no dofs it effectively has degree of 0
            IF (Face % BDOFs <= 0) CYCLE
            faceP = MAX(faceP, Face % PDefs % P)
         END DO
      END IF
!------------------------------------------------------------------------------
    END FUNCTION getFaceP
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
!> Get the number of Gauss points for the faces of p-elements
!------------------------------------------------------------------------------
    FUNCTION getNumberOfGaussPointsFace( Face, Mesh ) RESULT(ngp)
!------------------------------------------------------------------------------
      IMPLICIT NONE

      TYPE(Element_t), POINTER :: Face
      TYPE(Mesh_t) :: Mesh
      INTEGER :: ngp
!------------------------------------------------------------------------------
      TYPE(Element_t), POINTER :: Edge
      INTEGER :: edgeP, faceP, i, maxp
!------------------------------------------------------------------------------
      ! Get the max p of edges contained in the face
      edgeP = 0
      DO i=1, Face % TYPE % NumberOfEdges
         Edge => Mesh % Edges ( Face % EdgeIndexes(i) )
         edgeP = MAX(edgeP, Edge % PDefs % P)
      END DO

      IF (Face % BDOFs <= 0) THEN
        ! If no face dofs, the max p is defined by the edges
        maxp = edgeP
      ELSE
        maxp = MAX(edgeP, Face % PDefs % P)  
      END IF

      ! An economic quadrature may be available: 
      IF (Face % Pdefs % Serendipity .AND. Face % TYPE % ElementCode / 100 == 4) THEN
        !IF ( .NOT.(maxp < 4 .AND. Face % BDOFs>0) ) THEN
        IF (maxp > 1 .AND. maxp <= 8) THEN
          SELECT CASE(maxp)
          CASE(2)
            ngp = 8
          CASE(3)
            ngp = 12
          CASE(4)
            ngp = 20
          CASE(5)
            ngp = 25
          CASE(6)
            ngp = 36
          CASE(7)
            ngp = 45
          CASE(8)
            ngp = 60
          END SELECT
          RETURN
        END IF
      END IF

      ! Get the standard number of Gauss points:
      i = maxp + 1
      ngp = i ** 2
!------------------------------------------------------------------------------
    END FUNCTION getNumberOfGaussPointsFace
!------------------------------------------------------------------------------
 
  
!------------------------------------------------------------------------------
!>     Subroutine for getting reference p element nodes (because these are NOT
!>     yet defined in element description files)
!------------------------------------------------------------------------------
SUBROUTINE GetRefPElementNodes(Element, U, V, W)
        !------------------------------------------------------------------------------
        IMPLICIT NONE
        TYPE(ElementType_t) :: Element
        REAL(KIND=dp) :: U(:), V(:), W(:)
        !--------------------------------------------------------------------------------
        INTEGER :: n
        !--------------------------------------------------------------------------------    
        ! Reserve space for element nodes
        n = Element % NumberOfNodes

!$OMP single
        IF(.NOT.ALLOCATED(element % N_NodeU).AND.ALLOCATED(Element % NodeU) ) THEN
                ALLOCATE(Element % N_NodeU(n), Element % N_NodeV(n), Element % N_NodeW(n))
                element % N_NodeU = element % NodeU
                element % N_NodeV = element % NodeV
                element % N_NodeW = element % NOdeW
        END IF
!$OMP END single

    ! Select by element type given
    SELECT CASE(Element % ElementCode / 100)
    ! Line
    CASE(2)
       U(1:2) = (/ -1d0,1d0 /)
    ! Triangle
    CASE(3)
       U(1:3) = (/ -1d0,1d0,0d0 /)
       V(1:3) = (/ 0d0,0d0,SQRT(3.0d0) /)
    ! Quad
    CASE(4)
       U(1:4) = (/ -1d0,1d0,1d0,-1d0 /)
       V(1:4) = (/ -1d0,-1d0,1d0,1d0 /)
    ! Tetrahedron
    CASE(5)
       U(1:4) = (/ -1d0,1d0,0d0,0d0 /)
       V(1:4) = (/ 0d0,0d0,SQRT(3.0d0),1.0d0/SQRT(3.0d0) /)
       W(1:4) = (/ 0d0,0d0,0d0,2*SQRT(2.0d0/3.0d0) /)
    ! Pyramid
    CASE(6)
       U(1:5) = (/ -1d0,1d0,1d0,-1d0,0d0 /)
       V(1:5) = (/ -1d0,-1d0,1d0,1d0,0d0 /)
       W(1:5) = (/ 0d0,0d0,0d0,0d0,SQRT(2.0d0) /)
    ! Wedge
    CASE(7)
       U(1:6) = (/ -1d0,1d0,0d0,-1d0,1d0,0d0 /)
       V(1:6) = (/ 0d0,0d0,SQRT(3.0d0),0d0,0d0,SQRT(3.0d0) /)
       W(1:6) = (/ -1d0,-1d0,-1d0,1d0,1d0,1d0 /)
    ! Brick
    CASE(8)
       U(1:8) = (/ -1d0,1d0,1d0,-1d0,-1d0,1d0,1d0,-1d0 /)
       V(1:8) = (/ -1d0,-1d0,1d0,1d0,-1d0,-1d0,1d0,1d0 /)
       W(1:8) = (/ -1d0,-1d0,-1d0,-1d0,1d0,1d0,1d0,1d0 /)
    CASE DEFAULT
      WRITE(Message,'(A,I0)') 'Unknown element type: ',Element % ElementCode
      CALL Warn('PElementMaps::GetRefPElementNodes',Message)
    END SELECT
!------------------------------------------------------------------------------
    END SUBROUTINE GetRefPElementNodes
!------------------------------------------------------------------------------
    

END MODULE PElementMaps

!> \}
