diff --git a/src/modules/GlobalData/src/GlobalData.F90 b/src/modules/GlobalData/src/GlobalData.F90
index 522375441..9a6f4117d 100755
--- a/src/modules/GlobalData/src/GlobalData.F90
+++ b/src/modules/GlobalData/src/GlobalData.F90
@@ -388,7 +388,7 @@ MODULE GlobalData
 !! order = 2
 INTEGER(I4B), PARAMETER :: Quadrangle8 = 16
 !! order = 2, Serendipity
-INTEGER(I4B), PARAMETER :: Quadrangle16 = 160
+INTEGER(I4B), PARAMETER :: Quadrangle16 = 36
 !! order = 3
 !! Forder order 4, 5, and so on, we use 1601, 1602, 1603, etc.
 INTEGER(I4B), PARAMETER :: Tetrahedron = 4
diff --git a/src/submodules/Geometry/src/ReferenceElement_Method@ElementNameMethods.F90 b/src/submodules/Geometry/src/ReferenceElement_Method@ElementNameMethods.F90
index 48b637827..c4fec2589 100644
--- a/src/submodules/Geometry/src/ReferenceElement_Method@ElementNameMethods.F90
+++ b/src/submodules/Geometry/src/ReferenceElement_Method@ElementNameMethods.F90
@@ -142,7 +142,6 @@
 
 CASE DEFAULT
   ans = 0
-
 END SELECT
 END PROCEDURE Element_Type
 
@@ -193,6 +192,7 @@
   ans = ElementOrder_Pyramid(elemType)
 
 CASE DEFAULT
+
   ans = 0
 
 END SELECT
@@ -239,7 +239,8 @@
       Triangle21, &
       Quadrangle4, &
       Quadrangle8, &
-      Quadrangle9)
+      Quadrangle9, &
+      Quadrangle16)
   ans = 2
 CASE (Line2, &
       Line3, &
@@ -277,7 +278,7 @@
   & Triangle15b, Triangle21)
   ans = Triangle
 CASE (Quadrangle4, Quadrangle8, &
-  & Quadrangle9)
+  & Quadrangle9, Quadrangle16)
   ans = Quadrangle
 CASE (Tetrahedron4, Tetrahedron10, &
   & Tetrahedron20, Tetrahedron35, Tetrahedron56)
diff --git a/src/submodules/Geometry/src/ReferenceElement_Method@EnquireMethods.F90 b/src/submodules/Geometry/src/ReferenceElement_Method@EnquireMethods.F90
index b25a0d64b..6f9165be3 100644
--- a/src/submodules/Geometry/src/ReferenceElement_Method@EnquireMethods.F90
+++ b/src/submodules/Geometry/src/ReferenceElement_Method@EnquireMethods.F90
@@ -76,7 +76,8 @@
       Triangle21, &
       Quadrangle4, &
       Quadrangle8, &
-      Quadrangle9)
+      Quadrangle9, &
+      Quadrangle16)
   ans = .TRUE.
 CASE DEFAULT
   ans = .FALSE.
@@ -166,7 +167,7 @@
 
 MODULE PROCEDURE isQuadrangle1
 SELECT CASE (elemType)
-CASE (Quadrangle4, Quadrangle8, Quadrangle9)
+CASE (Quadrangle4, Quadrangle8, Quadrangle9, Quadrangle16)
   ans = .TRUE.
 CASE DEFAULT
   ans = .FALSE.
diff --git a/src/submodules/Geometry/src/ReferenceElement_Method@LocalNodeCoordsMethods.F90 b/src/submodules/Geometry/src/ReferenceElement_Method@LocalNodeCoordsMethods.F90
index a8efb7cf5..7475d8bf3 100644
--- a/src/submodules/Geometry/src/ReferenceElement_Method@LocalNodeCoordsMethods.F90
+++ b/src/submodules/Geometry/src/ReferenceElement_Method@LocalNodeCoordsMethods.F90
@@ -50,27 +50,27 @@
   ALLOCATE (nodecoord(3, 4))
   nodecoord = 0.0_DFP
   nodecoord(1, :) = [ &
-    & -1.0_DFP, 1.0_DFP, &
-    & -0.333333333333333_DFP, &
-    &  0.333333333333333_DFP]
+     -1.0_DFP, 1.0_DFP, &
+     -0.333333333333333_DFP, &
+      0.333333333333333_DFP]
 
 CASE (Line5)
   ALLOCATE (nodecoord(3, 5))
   nodecoord = 0.0_DFP
-  nodecoord(1, :) = [&
-    & -1.0_DFP, 1.0_DFP, &
-    & -0.5_DFP, 0.0_DFP, &
-    & 0.5_DFP]
+  nodecoord(1, :) = [ &
+     -1.0_DFP, 1.0_DFP, &
+     -0.5_DFP, 0.0_DFP, &
+     0.5_DFP]
 
 CASE (Line6)
   ALLOCATE (nodecoord(3, 6))
   nodecoord = 0.0_DFP
   nodecoord(1, :) = [ &
-    & -1.0_DFP, 1.0_DFP, &
-    & -0.666666666666666_DFP, &
-    & -0.333333333333333_DFP, &
-    & 0.666666666666666_DFP, &
-    & 0.333333333333333_DFP]
+     -1.0_DFP, 1.0_DFP, &
+     -0.666666666666666_DFP, &
+     -0.333333333333333_DFP, &
+     0.666666666666666_DFP, &
+     0.333333333333333_DFP]
 
 CASE (Triangle3)
   ALLOCATE (nodecoord(3, 3))
@@ -90,85 +90,85 @@
   ALLOCATE (nodecoord(3, 9))
   nodecoord = 0.0_DFP
   nodecoord(1, :) = [ &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.33333333333333333333_DFP, &
-    & 0.66666666666666666667_DFP, &
-    & 0.66666666666666666667_DFP, &
-    & 0.33333333333333333333_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP]
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.33333333333333333333_DFP, &
+     0.66666666666666666667_DFP, &
+     0.66666666666666666667_DFP, &
+     0.33333333333333333333_DFP, &
+     0.0_DFP, &
+     0.0_DFP]
 
   nodecoord(2, :) = [ &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.33333333333333333333_DFP, &
-    & 0.66666666666666666667_DFP, &
-    & 0.66666666666666666667_DFP, &
-    & 0.33333333333333333333_DFP]
+     0.0_DFP, &
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.33333333333333333333_DFP, &
+     0.66666666666666666667_DFP, &
+     0.66666666666666666667_DFP, &
+     0.33333333333333333333_DFP]
 
 CASE (Triangle10)
   ALLOCATE (nodecoord(3, 10))
   nodecoord = 0.0_DFP
   nodecoord(1, :) = [ &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.33333333333333333333_DFP, &
-    & 0.66666666666666666667_DFP, &
-    & 0.66666666666666666667_DFP, &
-    & 0.33333333333333333333_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.33333333333333333333_DFP]
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.33333333333333333333_DFP, &
+     0.66666666666666666667_DFP, &
+     0.66666666666666666667_DFP, &
+     0.33333333333333333333_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.33333333333333333333_DFP]
 
   nodecoord(2, :) = [ &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.33333333333333333333_DFP, &
-    & 0.66666666666666666667_DFP, &
-    & 0.66666666666666666667_DFP, &
-    & 0.33333333333333333333_DFP, &
-    & 0.33333333333333333333_DFP]
+     0.0_DFP, &
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.33333333333333333333_DFP, &
+     0.66666666666666666667_DFP, &
+     0.66666666666666666667_DFP, &
+     0.33333333333333333333_DFP, &
+     0.33333333333333333333_DFP]
 
 CASE (Triangle12)
   ! incomplete triangle; all nodes on boundary
   ALLOCATE (nodecoord(3, 12))
   nodecoord = 0.0_DFP
   nodecoord(1, :) = [ &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.25_DFP, &
-    & 0.5_DFP, &
-    & 0.75_DFP, &
-    & 0.75_DFP, &
-    & 0.5_DFP, &
-    & 0.25_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP]
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.25_DFP, &
+     0.5_DFP, &
+     0.75_DFP, &
+     0.75_DFP, &
+     0.5_DFP, &
+     0.25_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP]
 
   nodecoord(2, :) = [ &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.25_DFP, &
-    & 0.5_DFP, &
-    & 0.75_DFP, &
-    & 0.75_DFP, &
-    & 0.5_DFP, &
-    & 0.25_DFP]
+     0.0_DFP, &
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.25_DFP, &
+     0.5_DFP, &
+     0.75_DFP, &
+     0.75_DFP, &
+     0.5_DFP, &
+     0.25_DFP]
 
 CASE (Triangle15a)
   ! complete triangle; 12 nodes on boundary and
@@ -176,38 +176,38 @@
   ALLOCATE (nodecoord(3, 15))
   nodecoord = 0.0_DFP
   nodecoord(1, :) = [ &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.25_DFP, &
-    & 0.5_DFP, &
-    & 0.75_DFP, &
-    & 0.75_DFP, &
-    & 0.5_DFP, &
-    & 0.25_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.25_DFP, &
-    & 0.5_DFP, &
-    & 0.25_DFP]
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.25_DFP, &
+     0.5_DFP, &
+     0.75_DFP, &
+     0.75_DFP, &
+     0.5_DFP, &
+     0.25_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.25_DFP, &
+     0.5_DFP, &
+     0.25_DFP]
 
   nodecoord(2, :) = [ &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.25_DFP, &
-    & 0.5_DFP, &
-    & 0.75_DFP, &
-    & 0.75_DFP, &
-    & 0.5_DFP, &
-    & 0.25_DFP, &
-    & 0.25_DFP, &
-    & 0.25_DFP, &
-    & 0.5_DFP]
+     0.0_DFP, &
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.25_DFP, &
+     0.5_DFP, &
+     0.75_DFP, &
+     0.75_DFP, &
+     0.5_DFP, &
+     0.25_DFP, &
+     0.25_DFP, &
+     0.25_DFP, &
+     0.5_DFP]
 
 CASE (Triangle15b)
   ! Incomplete triangle
@@ -215,88 +215,88 @@
   nodecoord = 0.0_DFP
 
   nodecoord(1, :) = [ &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.2_DFP, &
-    & 0.4_DFP, &
-    & 0.6_DFP, &
-    & 0.8_DFP, &
-    & 0.8_DFP, &
-    & 0.6_DFP, &
-    & 0.4_DFP, &
-    & 0.2_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP]
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.2_DFP, &
+     0.4_DFP, &
+     0.6_DFP, &
+     0.8_DFP, &
+     0.8_DFP, &
+     0.6_DFP, &
+     0.4_DFP, &
+     0.2_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP]
 
   nodecoord(2, :) = [ &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.2_DFP, &
-    & 0.4_DFP, &
-    & 0.6_DFP, &
-    & 0.8_DFP, &
-    & 0.8_DFP, &
-    & 0.6_DFP, &
-    & 0.4_DFP, &
-    & 0.2_DFP]
+     0.0_DFP, &
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.2_DFP, &
+     0.4_DFP, &
+     0.6_DFP, &
+     0.8_DFP, &
+     0.8_DFP, &
+     0.6_DFP, &
+     0.4_DFP, &
+     0.2_DFP]
 
 CASE (Triangle21)
   ALLOCATE (nodecoord(3, 21))
   nodecoord = 0.0_DFP
 
   nodecoord(1, :) = [ &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.2_DFP, &
-    & 0.4_DFP, &
-    & 0.6_DFP, &
-    & 0.8_DFP, &
-    & 0.8_DFP, &
-    & 0.6_DFP, &
-    & 0.4_DFP, &
-    & 0.2_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.2_DFP, &
-    & 0.6_DFP, &
-    & 0.2_DFP, &
-    & 0.4_DFP, &
-    & 0.4_DFP, &
-    & 0.2_DFP]
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.2_DFP, &
+     0.4_DFP, &
+     0.6_DFP, &
+     0.8_DFP, &
+     0.8_DFP, &
+     0.6_DFP, &
+     0.4_DFP, &
+     0.2_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.2_DFP, &
+     0.6_DFP, &
+     0.2_DFP, &
+     0.4_DFP, &
+     0.4_DFP, &
+     0.2_DFP]
 
   nodecoord(2, :) = [ &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 1.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.0_DFP, &
-    & 0.2_DFP, &
-    & 0.4_DFP, &
-    & 0.6_DFP, &
-    & 0.8_DFP, &
-    & 0.8_DFP, &
-    & 0.6_DFP, &
-    & 0.4_DFP, &
-    & 0.2_DFP, &
-    & 0.2_DFP, &
-    & 0.2_DFP, &
-    & 0.6_DFP, &
-    & 0.2_DFP, &
-    & 0.4_DFP, &
-    & 0.4_DFP]
+     0.0_DFP, &
+     0.0_DFP, &
+     1.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.0_DFP, &
+     0.2_DFP, &
+     0.4_DFP, &
+     0.6_DFP, &
+     0.8_DFP, &
+     0.8_DFP, &
+     0.6_DFP, &
+     0.4_DFP, &
+     0.2_DFP, &
+     0.2_DFP, &
+     0.2_DFP, &
+     0.6_DFP, &
+     0.2_DFP, &
+     0.4_DFP, &
+     0.4_DFP]
 
 CASE (Quadrangle4)
   ALLOCATE (nodecoord(3, 4))
@@ -306,90 +306,110 @@
 
 CASE (Quadrangle8)
   nodecoord = RESHAPE([ &
-    & -1.0_DFP, -1.0_DFP, 0.0_DFP, &
-    & 1.0_DFP, -1.0_DFP, 0.0_DFP, &
-    & 1.0_DFP, 1.0_DFP, 0.0_DFP, &
-    & -1.0_DFP, 1.0_DFP, 0.0_DFP, &
-    & 0.0_DFP, -1.0_DFP, 0.0_DFP, &
-    & 1.0_DFP, 0.0_DFP, 0.0_DFP, &
-    & 0.0_DFP, 1.0_DFP, 0.0_DFP, &
-    & -1.0_DFP, 0.0_DFP, 0.0_DFP], [3, 8])
+     -1.0_DFP, -1.0_DFP, 0.0_DFP, &
+     1.0_DFP, -1.0_DFP, 0.0_DFP, &
+     1.0_DFP, 1.0_DFP, 0.0_DFP, &
+     -1.0_DFP, 1.0_DFP, 0.0_DFP, &
+     0.0_DFP, -1.0_DFP, 0.0_DFP, &
+     1.0_DFP, 0.0_DFP, 0.0_DFP, &
+     0.0_DFP, 1.0_DFP, 0.0_DFP, &
+     -1.0_DFP, 0.0_DFP, 0.0_DFP], [3, 8])
 
 CASE (Quadrangle9)
   nodecoord = RESHAPE([ &
-    & -1.0_DFP, -1.0_DFP, 0.0_DFP, &
-    & 1.0_DFP, -1.0_DFP, 0.0_DFP, &
-    & 1.0_DFP, 1.0_DFP, 0.0_DFP, &
-    & -1.0_DFP, 1.0_DFP, 0.0_DFP, &
-    & 0.0_DFP, -1.0_DFP, 0.0_DFP, &
-    & 1.0_DFP, 0.0_DFP, 0.0_DFP, &
-    & 0.0_DFP, 1.0_DFP, 0.0_DFP, &
-    & -1.0_DFP, 0.0_DFP, 0.0_DFP, &
-    & 0.0_DFP, 0.0_DFP, 0.0_DFP], [3, 9])
+     -1.0_DFP, -1.0_DFP, 0.0_DFP, &
+     1.0_DFP, -1.0_DFP, 0.0_DFP, &
+     1.0_DFP, 1.0_DFP, 0.0_DFP, &
+     -1.0_DFP, 1.0_DFP, 0.0_DFP, &
+     0.0_DFP, -1.0_DFP, 0.0_DFP, &
+     1.0_DFP, 0.0_DFP, 0.0_DFP, &
+     0.0_DFP, 1.0_DFP, 0.0_DFP, &
+     -1.0_DFP, 0.0_DFP, 0.0_DFP, &
+     0.0_DFP, 0.0_DFP, 0.0_DFP], [3, 9])
+
+CASE (Quadrangle16)
+  nodecoord = RESHAPE([ &
+     -1.0_DFP, -1.0_DFP, 0.0_DFP, &
+     1.0_DFP, -1.0_DFP, 0.0_DFP, &
+     1.0_DFP, 1.0_DFP, 0.0_DFP, &
+     -1.0_DFP, 1.0_DFP, 0.0_DFP, &
+     -1.0_dfp/3.0_dfp, -1.0_DFP, 0.0_DFP, &
+     1.0_dfp/3.0_dfp, -1.0_DFP, 0.0_DFP, &
+     1.0_dfp, -1.0_DFP/3.0_dfp, 0.0_DFP, &
+     1.0_dfp, 1.0_DFP/3.0_dfp, 0.0_DFP, &
+     1.0_dfp/3.0_dfp, 1.0_DFP, 0.0_DFP, &
+     -1.0_dfp/3.0_dfp, 1.0_DFP, 0.0_DFP, &
+     -1.0_dfp, 1.0_DFP/3.0_dfp, 0.0_DFP, &
+     -1.0_dfp, -1.0_DFP/3.0_dfp, 0.0_DFP, &
+     -1.0_dfp/3.0_Dfp, -1.0_DFP/3.0_dfp, 0.0_DFP, &
+     1.0_dfp/3.0_Dfp, -1.0_DFP/3.0_dfp, 0.0_DFP, &
+     1.0_dfp/3.0_Dfp, 1.0_DFP/3.0_dfp, 0.0_DFP, &
+     -1.0_dfp/3.0_Dfp, 1.0_DFP/3.0_dfp, 0.0_DFP], &
+    [3, 16])
 
 CASE (Hexahedron8)
   nodecoord = RESHAPE([ &
-  & -1.0_DFP, -1.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, -1.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, 1.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, 1.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, -1.0_DFP, 1.0_DFP, &
-  & 1.0_DFP, -1.0_DFP, 1.0_DFP, &
-  & 1.0_DFP, 1.0_DFP, 1.0_DFP, &
-  & -1.0_DFP, 1.0_DFP, 1.0_DFP], [3, 8])
+   -1.0_DFP, -1.0_DFP, -1.0_DFP, &
+   1.0_DFP, -1.0_DFP, -1.0_DFP, &
+   1.0_DFP, 1.0_DFP, -1.0_DFP, &
+   -1.0_DFP, 1.0_DFP, -1.0_DFP, &
+   -1.0_DFP, -1.0_DFP, 1.0_DFP, &
+   1.0_DFP, -1.0_DFP, 1.0_DFP, &
+   1.0_DFP, 1.0_DFP, 1.0_DFP, &
+   -1.0_DFP, 1.0_DFP, 1.0_DFP], [3, 8])
 
 CASE (Hexahedron20)
   nodecoord = RESHAPE([ &
-  & -1.0_DFP, -1.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, -1.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, 1.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, 1.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, -1.0_DFP, 1.0_DFP, &
-  & 1.0_DFP, -1.0_DFP, 1.0_DFP, &
-  & 1.0_DFP, 1.0_DFP, 1.0_DFP, &
-  & -1.0_DFP, 1.0_DFP, 1.0_DFP, &
-  & 0.0_DFP, -1.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, 0.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, -1.0_DFP, 0.0_DFP, &
-  & 1.0_DFP, 0.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, -1.0_DFP, 0.0_DFP, &
-  & 0.0_DFP, 1.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, 1.0_DFP, 0.0_DFP, &
-  & -1.0_DFP, 1.0_DFP, 0.0_DFP, &
-  & 0.0_DFP, -1.0_DFP, 1.0_DFP, &
-  & -1.0_DFP, 0.0_DFP, 1.0_DFP, &
-  & 1.0_DFP, 0.0_DFP, 1.0_DFP, &
-  & 0.0_DFP, 1.0_DFP, 1.0_DFP], [3, 20])
+   -1.0_DFP, -1.0_DFP, -1.0_DFP, &
+   1.0_DFP, -1.0_DFP, -1.0_DFP, &
+   1.0_DFP, 1.0_DFP, -1.0_DFP, &
+   -1.0_DFP, 1.0_DFP, -1.0_DFP, &
+   -1.0_DFP, -1.0_DFP, 1.0_DFP, &
+   1.0_DFP, -1.0_DFP, 1.0_DFP, &
+   1.0_DFP, 1.0_DFP, 1.0_DFP, &
+   -1.0_DFP, 1.0_DFP, 1.0_DFP, &
+   0.0_DFP, -1.0_DFP, -1.0_DFP, &
+   -1.0_DFP, 0.0_DFP, -1.0_DFP, &
+   -1.0_DFP, -1.0_DFP, 0.0_DFP, &
+   1.0_DFP, 0.0_DFP, -1.0_DFP, &
+   1.0_DFP, -1.0_DFP, 0.0_DFP, &
+   0.0_DFP, 1.0_DFP, -1.0_DFP, &
+   1.0_DFP, 1.0_DFP, 0.0_DFP, &
+   -1.0_DFP, 1.0_DFP, 0.0_DFP, &
+   0.0_DFP, -1.0_DFP, 1.0_DFP, &
+   -1.0_DFP, 0.0_DFP, 1.0_DFP, &
+   1.0_DFP, 0.0_DFP, 1.0_DFP, &
+   0.0_DFP, 1.0_DFP, 1.0_DFP], [3, 20])
 
 CASE (Hexahedron27)
   nodecoord = RESHAPE([ &
-  & -1.0_DFP, -1.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, -1.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, 1.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, 1.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, -1.0_DFP, 1.0_DFP, &
-  & 1.0_DFP, -1.0_DFP, 1.0_DFP, &
-  & 1.0_DFP, 1.0_DFP, 1.0_DFP, &
-  & -1.0_DFP, 1.0_DFP, 1.0_DFP, &
-  & 0.0_DFP, -1.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, 0.0_DFP, -1.0_DFP, &
-  & -1.0_DFP, -1.0_DFP, 0.0_DFP, &
-  & 1.0_DFP, 0.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, -1.0_DFP, 0.0_DFP, &
-  & 0.0_DFP, 1.0_DFP, -1.0_DFP, &
-  & 1.0_DFP, 1.0_DFP, 0.0_DFP, &
-  & -1.0_DFP, 1.0_DFP, 0.0_DFP, &
-  & 0.0_DFP, -1.0_DFP, 1.0_DFP, &
-  & -1.0_DFP, 0.0_DFP, 1.0_DFP, &
-  & 1.0_DFP, 0.0_DFP, 1.0_DFP, &
-  & 0.0_DFP, 1.0_DFP, 1.0_DFP, &
-  & 0.0_DFP, 0.0_DFP, -1.0_DFP, &
-  & 0.0_DFP, -1.0_DFP, 0.0_DFP, &
-  & -1.0_DFP, 0.0_DFP, 0.0_DFP, &
-  & 1.0_DFP, 0.0_DFP, 0.0_DFP, &
-  & 0.0_DFP, 1.0_DFP, 0.0_DFP, &
-  & 0.0_DFP, 0.0_DFP, 1.0_DFP, &
-  & 0.0_DFP, 0.0_DFP, 0.0_DFP], [3, 27])
+   -1.0_DFP, -1.0_DFP, -1.0_DFP, &
+   1.0_DFP, -1.0_DFP, -1.0_DFP, &
+   1.0_DFP, 1.0_DFP, -1.0_DFP, &
+   -1.0_DFP, 1.0_DFP, -1.0_DFP, &
+   -1.0_DFP, -1.0_DFP, 1.0_DFP, &
+   1.0_DFP, -1.0_DFP, 1.0_DFP, &
+   1.0_DFP, 1.0_DFP, 1.0_DFP, &
+   -1.0_DFP, 1.0_DFP, 1.0_DFP, &
+   0.0_DFP, -1.0_DFP, -1.0_DFP, &
+   -1.0_DFP, 0.0_DFP, -1.0_DFP, &
+   -1.0_DFP, -1.0_DFP, 0.0_DFP, &
+   1.0_DFP, 0.0_DFP, -1.0_DFP, &
+   1.0_DFP, -1.0_DFP, 0.0_DFP, &
+   0.0_DFP, 1.0_DFP, -1.0_DFP, &
+   1.0_DFP, 1.0_DFP, 0.0_DFP, &
+   -1.0_DFP, 1.0_DFP, 0.0_DFP, &
+   0.0_DFP, -1.0_DFP, 1.0_DFP, &
+   -1.0_DFP, 0.0_DFP, 1.0_DFP, &
+   1.0_DFP, 0.0_DFP, 1.0_DFP, &
+   0.0_DFP, 1.0_DFP, 1.0_DFP, &
+   0.0_DFP, 0.0_DFP, -1.0_DFP, &
+   0.0_DFP, -1.0_DFP, 0.0_DFP, &
+   -1.0_DFP, 0.0_DFP, 0.0_DFP, &
+   1.0_DFP, 0.0_DFP, 0.0_DFP, &
+   0.0_DFP, 1.0_DFP, 0.0_DFP, &
+   0.0_DFP, 0.0_DFP, 1.0_DFP, &
+   0.0_DFP, 0.0_DFP, 0.0_DFP], [3, 27])
 END SELECT
 
 END PROCEDURE Local_nodecoord
diff --git a/src/submodules/Geometry/src/ReferenceElement_Method@VTKMethods.F90 b/src/submodules/Geometry/src/ReferenceElement_Method@VTKMethods.F90
index e31a13535..17ecc9228 100644
--- a/src/submodules/Geometry/src/ReferenceElement_Method@VTKMethods.F90
+++ b/src/submodules/Geometry/src/ReferenceElement_Method@VTKMethods.F90
@@ -20,9 +20,10 @@
 ! summary: VTK format element
 
 SUBMODULE(ReferenceElement_Method) VTKMethods
+USE Utility, ONLY: arange
 IMPLICIT NONE
 
-INTEGER( I4B ), PARAMETER :: &
+INTEGER(I4B), PARAMETER :: &
   & vtk_point = 1, &
   & vtk_line2 = 3, &
   & vtk_triangle3 = 5, &
@@ -40,7 +41,10 @@
   & vtk_hexahedron27 = 29, &
   & vtk_prism15 = 26, &
   & vtk_prism18 = 32, &
-  & vtk_pyramid13 = 27
+  & vtk_line4 = 35, &
+  & vtk_pyramid13 = 27,  &
+  & vtk_quadrangle16 = 70
+! VTK_LAGRANGE_QUADRILATERAL
 CONTAINS
 
 !----------------------------------------------------------------------------
@@ -49,93 +53,102 @@
 
 MODULE PROCEDURE get_vtk_elemType
 
-SELECT CASE( ElemType )
-  CASE( Point1 )
-    vtk_type = vtk_point
-    nptrs = [1]
+SELECT CASE (ElemType)
+CASE (Point1)
+  vtk_type = vtk_point
+  nptrs = [1]
 
-  CASE( Line2 )
-    vtk_type = vtk_line2
-    nptrs = [1,2]
+CASE (Line2)
+  vtk_type = vtk_line2
+  nptrs = [1, 2]
 
-  CASE( Triangle3 )
-    vtk_type = vtk_triangle3
-    nptrs = [1,2,3]
+CASE (Triangle3)
+  vtk_type = vtk_triangle3
+  nptrs = [1, 2, 3]
 
-  CASE( Quadrangle4 )
-    vtk_type = vtk_quadrangle4
-    nptrs = [1,2,3,4]
+CASE (Quadrangle4)
+  vtk_type = vtk_quadrangle4
+  nptrs = [1, 2, 3, 4]
 
-  CASE( Tetrahedron4 )
-    vtk_type  = vtk_Tetrahedron4
-    nptrs = [1,2,3,4]
+CASE (Tetrahedron4)
+  vtk_type = vtk_Tetrahedron4
+  nptrs = [1, 2, 3, 4]
 
-  CASE( Hexahedron8 )
-    vtk_type  = vtk_Hexahedron8
-    nptrs = [1,2,3,4,5,6,7,8]
+CASE (Hexahedron8)
+  vtk_type = vtk_Hexahedron8
+  nptrs = [1, 2, 3, 4, 5, 6, 7, 8]
 
-  CASE( Prism6 )
-    vtk_type  = vtk_Prism6
-    nptrs = [1,2,3,4,5,6]
+CASE (Prism6)
+  vtk_type = vtk_Prism6
+  nptrs = [1, 2, 3, 4, 5, 6]
 
-  CASE( Pyramid5 )
-    vtk_type  = vtk_Pyramid5
-    nptrs = [1,2,3,4,5]
+CASE (Pyramid5)
+  vtk_type = vtk_Pyramid5
+  nptrs = [1, 2, 3, 4, 5]
 
   !! Order=2 elements
-  CASE( Line3 )
-    vtk_type  = vtk_line3
-    nptrs = [1,2,3]
-
-  CASE( Triangle6 )
-    vtk_type  = vtk_Triangle6
-    nptrs = [1,2,3,4,5,6]
-
-  CASE( Quadrangle9 )
-    vtk_type  = vtk_Quadrangle9
-    nptrs = [1,2,3,4,5,6,7,8,9]
-
-  CASE( Quadrangle8 )
-    vtk_type  = vtk_Quadrangle8
-    nptrs = [1,2,3,4,5,6,7,8]
-
-  CASE( Tetrahedron10 )
-    vtk_type  = vtk_Tetrahedron10
-    nptrs = 1+[0,1,2,3,4,5,6,7,9,8]
-
-  CASE( Hexahedron20 )
-    vtk_type  = vtk_Hexahedron20
-    nptrs = 1+[0,1,2,3,4,5,6,7, &
-      & 8,11,16,9,17,10,18,19,12,15,13,14]
-
-  CASE( Hexahedron27 )
-    vtk_type  = vtk_Hexahedron27
-    nptrs = 1+[0,1,2,3,4,5,6,7, &
-      & 8,11,16,9,17,10,18,19,12,15,13,14, &
-      & 24,22,20,21,23,25,26]
-
-  CASE( Prism15 )
-    vtk_type  = vtk_Prism15
-    nptrs = 1+[0,1,2,3,4,5, &
-      & 6,8,12,7,13,14,9,11,10]
-
-  CASE( Prism18 )
-    vtk_type  = vtk_Prism18
-    nptrs = 1+[0,1,2,3,4,5, &
-      & 6,8,12,7,13,14,9,11,10, &
-      & 15,17,16]
-
-  CASE( Pyramid13 )
-    vtk_type  = vtk_Pyramid13
-    nptrs = 1+[0,1,2,3,4,5, &
-      & 5,8,9,6,10,7,11,12]
-
-  CASE( Pyramid14 )
-  vtk_type  = vtk_Pyramid13
-    nptrs = 1+[0,1,2,3,4,5, &
-      & 5,8,9,6,10,7,11,12]
-
-  END SELECT
+CASE (Line3)
+  vtk_type = vtk_line3
+  nptrs = [1, 2, 3]
+
+CASE (Triangle6)
+  vtk_type = vtk_Triangle6
+  nptrs = [1, 2, 3, 4, 5, 6]
+
+CASE (Quadrangle9)
+  vtk_type = vtk_Quadrangle9
+  nptrs = [1, 2, 3, 4, 5, 6, 7, 8, 9]
+
+CASE (Quadrangle8)
+  vtk_type = vtk_Quadrangle8
+  nptrs = [1, 2, 3, 4, 5, 6, 7, 8]
+
+CASE (Tetrahedron10)
+  vtk_type = vtk_Tetrahedron10
+  nptrs = 1 + [0, 1, 2, 3, 4, 5, 6, 7, 9, 8]
+
+CASE (Hexahedron20)
+  vtk_type = vtk_Hexahedron20
+  nptrs = 1 + [0, 1, 2, 3, 4, 5, 6, 7, &
+    & 8, 11, 16, 9, 17, 10, 18, 19, 12, 15, 13, 14]
+
+CASE (Hexahedron27)
+  vtk_type = vtk_Hexahedron27
+  nptrs = 1 + [0, 1, 2, 3, 4, 5, 6, 7, &
+    & 8, 11, 16, 9, 17, 10, 18, 19, 12, 15, 13, 14, &
+    & 24, 22, 20, 21, 23, 25, 26]
+
+CASE (Prism15)
+  vtk_type = vtk_Prism15
+  nptrs = 1 + [0, 1, 2, 3, 4, 5, &
+    & 6, 8, 12, 7, 13, 14, 9, 11, 10]
+
+CASE (Prism18)
+  vtk_type = vtk_Prism18
+  nptrs = 1 + [0, 1, 2, 3, 4, 5, &
+    & 6, 8, 12, 7, 13, 14, 9, 11, 10, &
+    & 15, 17, 16]
+
+CASE (Pyramid13)
+  vtk_type = vtk_Pyramid13
+  nptrs = 1 + [0, 1, 2, 3, 4, 5, &
+    & 5, 8, 9, 6, 10, 7, 11, 12]
+
+CASE (Pyramid14)
+  vtk_type = vtk_Pyramid13
+  nptrs = 1 + [0, 1, 2, 3, 4, 5, &
+    & 5, 8, 9, 6, 10, 7, 11, 12]
+
+  !! order=3 element
+CASE (Line4)
+  vtk_type = vtk_line4
+  nptrs = [1, 2, 3, 4]
+
+CASE (Quadrangle16)
+  vtk_type = vtk_Quadrangle16
+  nptrs = [1, 2, 3, 4, 5, 6, 7, 8, 10, 9, &
+           12, 11, 13, 14, 16, 15]
+END SELECT
 END PROCEDURE get_vtk_elemType
 
-END SUBMODULE VTKMethods
\ No newline at end of file
+END SUBMODULE VTKMethods
diff --git a/src/submodules/Geometry/src/ReferenceQuadrangle_Method@Methods.F90 b/src/submodules/Geometry/src/ReferenceQuadrangle_Method@Methods.F90
index d5f95a52c..03b99b02d 100644
--- a/src/submodules/Geometry/src/ReferenceQuadrangle_Method@Methods.F90
+++ b/src/submodules/Geometry/src/ReferenceQuadrangle_Method@Methods.F90
@@ -585,8 +585,8 @@ END SUBROUTINE PARALLELOGRAMAREA2D
 DO iface = 1, 4
   DO ii = 1, order0 - 1
     con(2 + ii, iface) = jj + ii
-    jj = jj + 1
   END DO
+  jj = jj + order0 - 1
 END DO
 
 END PROCEDURE GetEdgeConnectivity_Quadrangle
diff --git a/src/submodules/Polynomial/src/QuadrangleInterpolationUtility@Methods.F90 b/src/submodules/Polynomial/src/QuadrangleInterpolationUtility@Methods.F90
index 10b233c25..f69357dcf 100644
--- a/src/submodules/Polynomial/src/QuadrangleInterpolationUtility@Methods.F90
+++ b/src/submodules/Polynomial/src/QuadrangleInterpolationUtility@Methods.F90
@@ -652,10 +652,7 @@
   ! internal nodes
   IF (ALL([p, q] .GE. 2_I4B)) THEN
 
-    CALL Reallocate( &
-      & xi_in,  &
-      & MAX(p - 1, 1_I4B), &
-      & MAX(q - 1_I4B, 1_I4B))
+    CALL Reallocate(xi_in, MAX(p - 1, 1_I4B), MAX(q - 1_I4B, 1_I4B))
     CALL Reallocate(eta_in, SIZE(xi_in, 1), SIZE(xi_in, 2))
     CALL Reallocate(temp_in, 2, SIZE(xi_in))
 
@@ -678,7 +675,7 @@
     xi_in = xi(ii1:ii2, jj1:jj2)
     eta_in = eta(ii1:ii2, jj1:jj2)
 
-    CALL IJ2VEFC_Quadrangle_Clockwise( &
+    CALL IJ2VEFC_Quadrangle_AntiClockwise( &
       & xi=xi_in,  &
       & eta=eta_in, &
       & temp=temp_in, &