diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90
index 7a77ddc..5a9e512 100644
--- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90
+++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90
@@ -4,9 +4,11 @@
 !              z == unused
 MODULE moduleMesh1DCart
   USE moduleMesh
-  use moduleMeshCommon
   IMPLICIT NONE
 
+  REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0,       DSQRT(3.D0/5.D0) /)
+  REAL(8), PARAMETER:: wSeg(1:3)   = (/        5.D0/9.D0 , 8.D0/9.D0,        5.D0/9.D0  /)
+
   TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DCart
     !Element coordinates
     REAL(8):: x = 0.D0
@@ -164,6 +166,18 @@ MODULE moduleMesh1DCart
 
     END FUNCTION randPosEdge
 
+    !Compute element functions at point Xi
+    pure function fPsiPoint(Xi, nNodes) RESULT(fPsi)
+      implicit none
+
+      real(8), intent(in):: Xi(1:3)
+      integer, intent(in):: nNodes
+      real(8):: fPsi(1:nNodes)
+
+      fPsi = (/ 1.D0 /)
+
+    end function fPsiPoint
+
     !CELL FUNCTIONS
     !SEGMENT FUNCTIONS
     !Init element
@@ -228,6 +242,35 @@ MODULE moduleMesh1DCart
 
     END FUNCTION randPos1DCartSegm
 
+    !Compute element functions at point Xi
+    PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi = (/ 1.D0 - Xi(1), &
+                1.D0 + Xi(1) /)
+
+      fPsi    = fPsi * 0.50D0
+
+    END FUNCTION fPsiSegm
+
+    !Derivative element function at coordinates Xi
+    PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: dPsi(1:3,1:nNodes)
+
+      dPsi = 0.D0
+
+      dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /)
+
+    END FUNCTION dPsiSegm
+
     !Partial derivative in global coordinates
     PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer)
       IMPLICIT NONE
diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90
index fcc438a..6e35fe2 100644
--- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90
+++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90
@@ -4,9 +4,11 @@
 !              z == unused
 MODULE moduleMesh1DRad
   USE moduleMesh
-  use moduleMeshCommon
   IMPLICIT NONE
 
+  REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0,       DSQRT(3.D0/5.D0) /)
+  REAL(8), PARAMETER:: wSeg(1:3)   = (/        5.D0/9.D0 , 8.D0/9.D0,        5.D0/9.D0  /)
+
   TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DRad
     !Element coordinates
     REAL(8):: r = 0.D0
@@ -164,6 +166,18 @@ MODULE moduleMesh1DRad
 
     END FUNCTION randPos1DRad
 
+    !Compute element functions at point Xi
+    pure function fPsiPoint(Xi, nNodes) RESULT(fPsi)
+      implicit none
+
+      real(8), intent(in):: Xi(1:3)
+      integer, intent(in):: nNodes
+      real(8):: fPsi(1:nNodes)
+
+      fPsi = (/ 1.D0 /)
+
+    end function fPsiPoint
+
     !CELL FUNCTIONS
     !SEGMENT FUNCTIONS
     !Init element
@@ -228,6 +242,35 @@ MODULE moduleMesh1DRad
 
     END FUNCTION randPos1DRadSegm
 
+    !Compute element functions at point Xi
+    PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi = (/ 1.D0 - Xi(1), &
+                1.D0 + Xi(1) /)
+
+      fPsi    = fPsi * 0.50D0
+
+    END FUNCTION fPsiSegm
+
+    !Derivative element function at coordinates Xi
+    PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: dPsi(1:3,1:nNodes)
+
+      dPsi = 0.D0
+
+      dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /)
+
+    END FUNCTION dPsiSegm
+
     !Partial derivative in global coordinates
     PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer)
       IMPLICIT NONE
diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
index 5bfb1e4..e2014a7 100644
--- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90
+++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
@@ -4,9 +4,16 @@
 !                  z == unused
 MODULE moduleMesh2DCart
   USE moduleMesh
-  use moduleMeshCommon
   IMPLICIT NONE
 
+  !Values for Gauss integral
+  REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0,      DSQRT(3.D0/5.D0) /)
+  REAL(8), PARAMETER:: wQuad(1:3)   = (/        5.D0/9.D0,  8.D0/9.D0,       5.D0/9.D0  /)
+
+  REAL(8), PARAMETER:: Xi1Tria(1:4) = (/   1.D0/3.D0,   1.D0/5.D0,   3.D0/5.D0,   1.D0/5.D0  /)
+  REAL(8), PARAMETER:: Xi2Tria(1:4) = (/   1.D0/3.D0,   1.D0/5.D0,   1.D0/5.D0,   3.D0/5.D0  /)
+  REAL(8), PARAMETER:: wTria(1:4)   = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /)
+
   TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCart
     !Element coordinates
     REAL(8):: x = 0.D0, y = 0.D0
@@ -221,6 +228,21 @@ MODULE moduleMesh2DCart
 
     END FUNCTION randPosEdge
 
+    !Compute element functions at point Xi
+    pure function fPsiSegm(Xi, nNodes) RESULT(fPsi)
+      implicit none
+
+      real(8), intent(in):: Xi(1:3)
+      integer, intent(in):: nNodes
+      real(8):: fPsi(1:nNodes)
+
+      fPsi = (/ 1.D0 - Xi(1), &
+                1.D0 + Xi(1) /)
+
+      fPsi    = fPsi * 0.50D0
+
+    end function fPsiSegm
+
     !VOLUME FUNCTIONS
     !QUAD FUNCTIONS
     !Init element
@@ -297,6 +319,49 @@ MODULE moduleMesh2DCart
 
     END FUNCTION randPosCellQuad
 
+    !Compute element functions in point Xi
+    PURE FUNCTION fPsiQuad(Xi, nNodes) RESULT(fPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi = 0.D0
+
+      fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
+                (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
+                (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
+                (1.D0 - Xi(1)) * (1.D0 + Xi(2)) /)
+
+      fPsi    = fPsi * 0.25D0
+
+    END FUNCTION fPsiQuad
+
+    !Derivative element function at coordinates Xi
+    PURE FUNCTION dPsiQuad(Xi, nNodes) RESULT(dPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: dPsi(1:3,1:nNodes)
+
+      dPsi = 0.D0
+
+      dPsi(1, 1:4) = (/ -(1.D0 - Xi(2)), &
+                         (1.D0 - Xi(2)), &
+                         (1.D0 + Xi(2)), &
+                        -(1.D0 + Xi(2)) /)
+
+      dPsi(2, 1:4) = (/ -(1.D0 - Xi(1)), &
+                        -(1.D0 + Xi(1)), &
+                         (1.D0 + Xi(1)), &
+                         (1.D0 - Xi(1)) /)
+
+      dPsi = dPsi * 0.25D0
+
+    END FUNCTION dPsiQuad
+
     !Partial derivative in global coordinates
     PURE FUNCTION partialDerQuad(self, nNodes, dPsi) RESULT(pDer)
       IMPLICIT NONE
@@ -600,6 +665,35 @@ MODULE moduleMesh2DCart
 
     END FUNCTION randPosCellTria
 
+    !Compute element functions in point Xi
+    PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi(1) = 1.D0 - Xi(1) - Xi(2)
+      fPsi(2) =        Xi(1)
+      fPsi(3) =                Xi(2)
+
+    END FUNCTION fPsiTria
+
+    !Compute element derivative functions in point Xi
+    PURE FUNCTION dPsiTria(Xi, nNodes) RESULT(dPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: dPsi(1:3,1:nNodes)
+
+      dPsi = 0.D0
+
+      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
+      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
+
+    END FUNCTION dPsiTria
+
     !Compute the derivatives in global coordinates
     PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer)
       IMPLICIT NONE
diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
index e606641..8b64b9d 100644
--- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
+++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
@@ -4,9 +4,16 @@
 !               z == theta (unused)
 MODULE moduleMesh2DCyl
   USE moduleMesh
-  use moduleMeshCommon
   IMPLICIT NONE
 
+  !Values for Gauss integral
+  REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0,      DSQRT(3.D0/5.D0) /)
+  REAL(8), PARAMETER:: wQuad(1:3)   = (/        5.D0/9.D0,  8.D0/9.D0,       5.D0/9.D0  /)
+
+  REAL(8), PARAMETER:: Xi1Tria(1:4) = (/   1.D0/3.D0,   1.D0/5.D0,   3.D0/5.D0,   1.D0/5.D0  /)
+  REAL(8), PARAMETER:: Xi2Tria(1:4) = (/   1.D0/3.D0,   1.D0/5.D0,   1.D0/5.D0,   3.D0/5.D0  /)
+  REAL(8), PARAMETER:: wTria(1:4)   = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /)
+
   TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCyl
     !Element coordinates
     REAL(8):: r = 0.D0, z = 0.D0
@@ -230,6 +237,21 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION randPosEdge
 
+    !Compute element functions at point Xi
+    pure function fPsiSegm(Xi, nNodes) RESULT(fPsi)
+      implicit none
+
+      real(8), intent(in):: Xi(1:3)
+      integer, intent(in):: nNodes
+      real(8):: fPsi(1:nNodes)
+
+      fPsi = (/ 1.D0 - Xi(1), &
+                1.D0 + Xi(1) /)
+
+      fPsi    = fPsi * 0.50D0
+
+    end function fPsiSegm
+
     !CELL FUNCTIONS
     !QUAD FUNCTIONS
     !Init element
@@ -305,6 +327,49 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION randPosCellQuad
 
+    !Computes element functions in point Xi
+    PURE FUNCTION fPsiQuad(Xi, nNodes) RESULT(fPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi = 0.D0
+
+      fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
+                (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
+                (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
+                (1.D0 - Xi(1)) * (1.D0 + Xi(2)) /)
+
+      fPsi    = fPsi * 0.25D0
+
+    END FUNCTION fPsiQuad
+
+    !Derivative element function at coordinates Xi
+    PURE FUNCTION dPsiQuad(Xi, nNodes) RESULT(dPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: dPsi(1:3,1:nNodes)
+
+      dPsi = 0.D0
+
+      dPsi(1, 1:4) = (/ -(1.D0 - Xi(2)), &
+                         (1.D0 - Xi(2)), &
+                         (1.D0 + Xi(2)), &
+                        -(1.D0 + Xi(2)) /)
+
+      dPsi(2, 1:4) = (/ -(1.D0 - Xi(1)), &
+                        -(1.D0 + Xi(1)), &
+                         (1.D0 + Xi(1)), &
+                         (1.D0 - Xi(1)) /)
+
+      dPsi = dPsi * 0.25D0
+
+    END FUNCTION dPsiQuad
+
     !Partial derivative in global coordinates
     PURE FUNCTION partialDerQuad(self, nNodes, dPsi) RESULT(pDer)
       IMPLICIT NONE
@@ -627,6 +692,35 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION randPosCellTria
 
+    !Compute element functions in point Xi
+    PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi(1) = 1.D0 - Xi(1) - Xi(2)
+      fPsi(2) =        Xi(1)
+      fPsi(3) =                Xi(2)
+
+    END FUNCTION fPsiTria
+
+    !Compute element derivative functions in point Xi
+    PURE FUNCTION dPsiTria(Xi, nNodes) RESULT(dPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: dPsi(1:3,1:nNodes)
+
+      dPsi = 0.D0
+
+      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
+      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
+
+    END FUNCTION dPsiTria
+
     !Compute the derivatives in global coordinates
     PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer)
       IMPLICIT NONE
diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90
index 1655ab8..a10616a 100644
--- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90
+++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90
@@ -4,7 +4,6 @@
 !                  z == z
 MODULE moduleMesh3DCart
   USE moduleMesh
-  use moduleMeshCommon
   IMPLICIT NONE
 
   TYPE, PUBLIC, EXTENDS(meshNode):: meshNode3DCart
@@ -205,6 +204,20 @@ MODULE moduleMesh3DCart
 
     END FUNCTION randPosEdgeTria
 
+    !Shape functions for triangular surface
+    PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi(1) = 1.D0 - Xi(1) - Xi(2)
+      fPsi(2) =        Xi(1)
+      fPsi(3) =                Xi(2)
+
+    END FUNCTION fPsiTria
+
     !VOLUME FUNCTIONS
     !TETRA FUNCTIONS
     !Init element
@@ -282,6 +295,37 @@ MODULE moduleMesh3DCart
 
     END FUNCTION randPosCellTetra
 
+    !Compute element functions in point Xi
+    PURE FUNCTION fPsiTetra(Xi, nNodes) RESULT(fPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi(1) = 1.D0 - Xi(1) - Xi(2) - Xi(3)
+      fPsi(2) =        Xi(1)
+      fPsi(3) =                Xi(2)
+      fPsi(4) =                        Xi(3)
+
+    END FUNCTION fPsiTetra
+
+    !Compute element derivative functions in point Xi
+    PURE FUNCTION dPsiTetra(Xi, nNodes) RESULT(dPsi)
+      IMPLICIT NONE
+
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: dPsi(1:3, 1:nNodes)
+
+      dPsi = 0.D0
+
+      dPsi(1,1:4) = (/ -1.D0, 1.D0, 0.D0, 0.D0 /)
+      dPsi(2,1:4) = (/ -1.D0, 0.D0, 1.D0, 0.D0 /)
+      dPsi(3,1:4) = (/ -1.D0, 0.D0, 0.D0, 1.D0 /)
+
+    END FUNCTION dPsiTetra
+
     !Compute the derivatives in global coordinates
     PURE FUNCTION partialDerTetra(self, nNodes, dPsi) RESULT(pDer)
       IMPLICIT NONE
diff --git a/src/modules/mesh/makefile b/src/modules/mesh/makefile
index 7481690..e21187f 100644
--- a/src/modules/mesh/makefile
+++ b/src/modules/mesh/makefile
@@ -18,7 +18,7 @@ all: moduleMesh.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o
 0D.o: moduleMesh.o
 	$(MAKE) -C 0D all
 
-moduleMesh.o: moduleMeshCommon.o moduleMesh.f90
+moduleMesh.o: moduleMesh.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 	$(FC) $(FCFLAGS) -c moduleMesh@elements.f90 -o $(OBJDIR)/$@
 	$(FC) $(FCFLAGS) -c moduleMesh@boundary.f90 -o $(OBJDIR)/$@
diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundary.f90
similarity index 99%
rename from src/modules/mesh/moduleMesh@boundaryParticle.f90
rename to src/modules/mesh/moduleMesh@boundary.f90
index a6ce1b1..cf7a45b 100644
--- a/src/modules/mesh/moduleMesh@boundaryParticle.f90
+++ b/src/modules/mesh/moduleMesh@boundary.f90
@@ -1,5 +1,5 @@
 !moduleMeshBoundary: Boundary functions for the mesh edges
-submodule(moduleMesh) boundaryParticle
+submodule(moduleMesh) boundary
   contains
     module function boundaryParticleName_to_Index(boundaryName) result(bp)
       use moduleErrors
@@ -502,4 +502,4 @@ submodule(moduleMesh) boundaryParticle
 
     end function physicalSurface_to_id
 
-end submodule boundaryParticle
+end submodule boundary
diff --git a/src/modules/mesh/moduleMesh@common.f90 b/src/modules/mesh/moduleMesh@common.f90
new file mode 100644
index 0000000..fe710a8
--- /dev/null
+++ b/src/modules/mesh/moduleMesh@common.f90
@@ -0,0 +1,3 @@
+module moduleMeshCommon
+
+end module moduleMeshCommon
diff --git a/src/modules/mesh/moduleMeshCommon.f90 b/src/modules/mesh/moduleMeshCommon.f90
deleted file mode 100644
index d2103b0..0000000
--- a/src/modules/mesh/moduleMeshCommon.f90
+++ /dev/null
@@ -1,164 +0,0 @@
-! Common subroutines and parameters for mesh elements
-module moduleMeshCommon
-  ! values for integration gauss integral
-  ! Segment
-  real(8), parameter:: corSeg(1:3)  = (/ -dsqrt(3.D0/5.D0),  0.D0,      dsqrt(3.D0/5.D0) /)
-  real(8), parameter:: wSeg(1:3)    = (/        5.D0/9.D0,   8.D0/9.D0,       5.D0/9.D0  /)
-
-  ! Quad
-  real(8), parameter:: corQuad(1:3) = (/ -dsqrt(3.D0/5.D0),  0.D0,      dsqrt(3.D0/5.D0) /)
-  real(8), parameter:: wQuad(1:3)   = (/        5.D0/9.D0,   8.D0/9.D0,       5.D0/9.D0  /)
-
-  ! tria
-  real(8), parameter:: Xi1Tria(1:4) = (/        1.D0/3.D0,   1.D0/5.D0,       3.D0/5.D0,   1.D0/5.D0  /)
-  real(8), parameter:: Xi2Tria(1:4) = (/        1.D0/3.D0,   1.D0/5.D0,       1.D0/5.D0,   3.D0/5.D0  /)
-  real(8), parameter:: wTria(1:4)   = (/      -27.D0/96.D0, 25.D0/96.D0,     25.D0/96.D0, 25.D0/96.D0 /)
-
-  contains
-    ! ELEMENT FUNCTIONS
-    ! Point
-    pure function fPsiPoint(Xi, nNodes) RESULT(fPsi)
-      implicit none
-
-      real(8), intent(in):: Xi(1:3)
-      integer, intent(in):: nNodes
-      real(8):: fPsi(1:nNodes)
-
-      fPsi = (/ 1.D0 /)
-
-    end function fPsiPoint
-
-    ! Segment
-    pure function fPsiSegm(Xi, nNodes) result(fPsi)
-      implicit none
-
-      real(8), intent(in):: Xi(1:3)
-      integer, intent(in):: nNodes
-      real(8):: fPsi(1:nNodes)
-
-      fPsi = (/ 1.D0 - Xi(1), &
-                1.D0 + Xi(1) /)
-
-      fPsi    = fPsi * 0.50D0
-
-    end function fPsiSegm
-
-    ! Quad
-    pure function fPsiQuad(Xi, nNodes) result(fPsi)
-      implicit none
-
-      real(8), intent(in):: Xi(1:3)
-      integer, intent(in):: nNodes
-      real(8):: fPsi(1:nNodes)
-
-      fPsi = 0.D0
-
-      fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
-                (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
-                (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
-                (1.D0 - Xi(1)) * (1.D0 + Xi(2)) /)
-
-      fPsi    = fPsi * 0.25D0
-
-    end function fPsiQuad
-
-    ! Tria
-    pure function fPsiTria(Xi, nNodes) result(fPsi)
-      IMPLICIT NONE
-
-      real(8), intent(in):: Xi(1:3)
-      integer, intent(in):: nNodes
-      real(8):: fPsi(1:nNodes)
-
-      fPsi(1) = 1.D0 - Xi(1) - Xi(2)
-      fPsi(2) =        Xi(1)
-      fPsi(3) =                Xi(2)
-
-    end function fPsiTria
-
-    ! Tetra
-    pure function fPsiTetra(Xi, nNodes) result(fPsi)
-      implicit none
-
-      real(8), intent(in):: Xi(1:3)
-      integer, intent(in):: nNodes
-      real(8):: fPsi(1:nNodes)
-
-      fPsi(1) = 1.D0 - Xi(1) - Xi(2) - Xi(3)
-      fPsi(2) =        Xi(1)
-      fPsi(3) =                Xi(2)
-      fPsi(4) =                        Xi(3)
-
-    end function fPsiTetra
-
-    ! DERIVATIVE FOR ELEMENT FUNCTIONS
-    ! Segment
-    pure function dPsiSegm(Xi, nNodes) result(dPsi)
-      implicit none
-
-      real(8), intent(in):: Xi(1:3)
-      integer, intent(in):: nNodes
-      real(8):: dPsi(1:3,1:nNodes)
-
-      dPsi = 0.D0
-
-      dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /)
-
-    end function dPsiSegm
-
-    ! Quad
-    pure function dPsiQuad(Xi, nNodes) result(dPsi)
-      implicit none
-
-      real(8), intent(in):: Xi(1:3)
-      integer, intent(in):: nNodes
-      real(8):: dPsi(1:3,1:nNodes)
-
-      dPsi = 0.D0
-
-      dPsi(1, 1:4) = (/ -(1.D0 - Xi(2)), &
-                         (1.D0 - Xi(2)), &
-                         (1.D0 + Xi(2)), &
-                        -(1.D0 + Xi(2)) /)
-
-      dPsi(2, 1:4) = (/ -(1.D0 - Xi(1)), &
-                        -(1.D0 + Xi(1)), &
-                         (1.D0 + Xi(1)), &
-                         (1.D0 - Xi(1)) /)
-
-      dPsi = dPsi * 0.25D0
-
-    end function dPsiQuad
-
-    ! Tria
-    pure function dPsiTria(Xi, nNodes) result(dPsi)
-      implicit none
-
-      real(8), intent(in):: Xi(1:3)
-      integer, intent(in):: nNodes
-      real(8):: dPsi(1:3,1:nNodes)
-
-      dPsi = 0.D0
-
-      dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /)
-      dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /)
-
-    end function dPsiTria
-
-    !Compute element derivative functions in point Xi
-    pure function dPsiTetra(Xi, nNodes) result(dPsi)
-      implicit none
-
-      real(8), intent(in):: Xi(1:3)
-      integer, intent(in):: nNodes
-      real(8):: dPsi(1:3, 1:nNodes)
-
-      dPsi = 0.D0
-
-      dPsi(1,1:4) = (/ -1.D0, 1.D0, 0.D0, 0.D0 /)
-      dPsi(2,1:4) = (/ -1.D0, 0.D0, 1.D0, 0.D0 /)
-      dPsi(3,1:4) = (/ -1.D0, 0.D0, 0.D0, 1.D0 /)
-
-    end function dPsiTetra
-
-end module moduleMeshCommon