DOES NOT COMPILE: Break

Small break of changing functions.
Still some geometries to change.

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -19,6 +19,7 @@ MODULE moduleMesh2DCyl
     !Element coordinates
     REAL(8):: r = 0.D0, z = 0.D0
     CONTAINS
+      !meshNode DEFERRED PROCEDURES
       PROCEDURE, PASS:: init           => initNode2DCyl
       PROCEDURE, PASS:: getCoordinates => getCoord2DCyl
 
@@ -30,6 +31,7 @@ MODULE moduleMesh2DCyl
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
     CONTAINS
+      !meshEdge DEFERRED PROCEDURES
       PROCEDURE, PASS:: init         => initEdge2DCyl
       PROCEDURE, PASS:: getNodes     => getNodes2DCyl
       PROCEDURE, PASS:: intersection => intersection2DCylEdge
@@ -37,28 +39,8 @@ MODULE moduleMesh2DCyl
 
   END TYPE meshEdge2DCyl
 
-  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell2DCyl
-    CONTAINS
-      PROCEDURE, PASS:: detJac    => detJ2DCyl
-      PROCEDURE, PASS:: invJac    => invJ2DCyl
-      PROCEDURE(partialDer_interface), DEFERRED, PASS, PRIVATE:: partialDer
-
-  END TYPE meshCell2DCyl
-
-  ABSTRACT INTERFACE
-    PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dz, dr)
-      IMPORT meshCell2DCyl
-      CLASS(meshCell2DCyl), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
-      REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
-
-    END SUBROUTINE partialDer_interface
-
-  END INTERFACE
-
   !Quadrilateral volume element
-  TYPE, PUBLIC, EXTENDS(meshCell2DCyl):: meshCell2DCylQuad
+  TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylQuad
     !Element coordinates
     REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0
     !Connectivity to nodes
@@ -68,25 +50,29 @@ MODULE moduleMesh2DCyl
     REAL(8):: arNodes(1:4) = 0.D0
 
     CONTAINS
+      !meshCell DEFERRED PROCEDURES
       PROCEDURE, PASS::   init                => initCellQuad2DCyl
+      PROCEDURE, PASS::   getNodes            => getNodesQuad
       PROCEDURE, PASS::   randPos             => randPosCellQuad
-      PROCEDURE, PASS::   area        => areaQuad
+      PROCEDURE, NOPASS:: fPsi                => fPsiQuad
-      PROCEDURE, PASS::   fPsi        => fPsiQuad
+      PROCEDURE, NOPASS:: dPsi                => dPsiQuad
-      PROCEDURE, PASS::   dPsi        => dPsiQuad
+      PROCEDURE, PASS::   partialDer          => partialDerQuad
-      PROCEDURE, PASS, PRIVATE::   partialDer  => partialDerQuad
+      PROCEDURE, NOPASS:: detJac              => detJ2DCyl
-      PROCEDURE, PASS::   elemK       => elemKQuad
+      PROCEDURE, NOPASS:: invJac              => invJ2DCyl
-      PROCEDURE, PASS::   elemF       => elemFQuad
       PROCEDURE, PASS::   gatherElectricField => gatherEFQuad
       PROCEDURE, PASS::   gatherMagneticField => gatherMFQuad
+      PROCEDURE, PASS::   elemK               => elemKQuad
+      PROCEDURE, PASS::   elemF               => elemFQuad
       PROCEDURE, NOPASS:: inside              => insideQuad
-      PROCEDURE, PASS::   getNodes    => getNodesQuad
       PROCEDURE, PASS::   phy2log             => phy2logQuad
-      PROCEDURE, PASS::   nextElement => nextElementQuad
+      PROCEDURE, PASS::   neighbourElement    => neighbourElementQuad
+      !PARTICLUAR PROCEDURES
+      PROCEDURE, PASS::   area                => areaQuad
 
   END TYPE meshCell2DCylQuad
 
   !Triangular volume element
-  TYPE, PUBLIC, EXTENDS(meshCell2DCyl):: meshCell2DCylTria
+  TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylTria
     !Element coordinates
     REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0
     !Connectivity to nodes
@@ -96,20 +82,24 @@ MODULE moduleMesh2DCyl
     REAL(8):: arNodes(1:3) = 0.D0
 
     CONTAINS
+      !meshCell DEFERRED PROCEDURES
       PROCEDURE, PASS::   init                => initCellTria2DCyl
+      PROCEDURE, PASS::   getNodes            => getNodesTria
       PROCEDURE, PASS::   randPos             => randPosCellTria
-      PROCEDURE, PASS::   area        => areaTria
+      PROCEDURE, NOPASS:: fPsi                => fPsiTria
-      PROCEDURE, PASS::   fPsi        => fPsiTria
+      PROCEDURE, NOPASS:: dPsi                => dPsiTria
-      PROCEDURE, PASS::   dPsi        => dPsiTria
+      PROCEDURE, PASS::   partialDer          => partialDerTria
-      PROCEDURE, PASS, PRIVATE::   partialDer  => partialDerTria
+      PROCEDURE, NOPASS:: detJac              => detJ2DCyl
-      PROCEDURE, PASS::   elemK       => elemKTria
+      PROCEDURE, NOPASS:: invJac              => invJ2DCyl
-      PROCEDURE, PASS::   elemF       => elemFTria
       PROCEDURE, PASS::   gatherElectricField => gatherEFTria
       PROCEDURE, PASS::   gatherMagneticField => gatherMFTria
+      PROCEDURE, PASS::   elemK               => elemKTria
+      PROCEDURE, PASS::   elemF               => elemFTria
       PROCEDURE, NOPASS:: inside              => insideTria
-      PROCEDURE, PASS::   getNodes    => getNodesTria
       PROCEDURE, PASS::   phy2log             => phy2logTria
-      PROCEDURE, PASS::   nextElement => nextElementTria
+      PROCEDURE, PASS::   neighbourElement    => neighbourElementTria
+      !PARTICULAR PROCEDURES
+      PROCEDURE, PASS::   area                => areaTria
 
   END TYPE meshCell2DCylTria
 
@@ -294,99 +284,17 @@ MODULE moduleMesh2DCyl
 
     END SUBROUTINE initCellQuad2DCyl
 
-    !Computes element area
+    !Gets nodes from quadrilateral element
-    PURE SUBROUTINE areaQuad(self)
+    PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n)
-      USE moduleConstParam, ONLY: PI8
-      IMPLICIT NONE
-
-      CLASS(meshCell2DCylQuad), INTENT(inout):: self
-      REAL(8):: r, Xi(1:3)
-      REAL(8):: detJ
-      REAL(8):: fPsi(1:4), fPsi_node(1:4)
-
-      self%volume  = 0.D0
-      self%arNodes = 0.D0
-      !2D 1 point Gauss Quad Integral
-      Xi = 0.D0
-      detJ = self%detJac(Xi, 4)*PI8 !4*2*pi
-      fPsi = self%fPsi(Xi, 4)
-      !Computes total volume of the cell
-      r = DOT_PRODUCT(fPsi,self%r)
-      self%volume  = r*detJ
-      !Computes volume per node
-      Xi = (/-5.D-1, -5.D-1, 0.D0/)
-      r = self%gatherF(Xi, 4, self%r)
-      self%arNodes(1) = fPsi(1)*r*detJ
-      Xi = (/ 5.D-1, -5.D-1, 0.D0/)
-      r = self%gatherF(Xi, 4, self%r)
-      self%arNodes(2) = fPsi(2)*r*detJ
-      Xi = (/ 5.D-1,  5.D-1, 0.D0/)
-      r = self%gatherF(Xi, 4, self%r)
-      self%arNodes(3) = fPsi(3)*r*detJ
-      Xi = (/-5.D-1,  5.D-1, 0.D0/)
-      r = self%gatherF(Xi, 4, self%r)
-      self%arNodes(4) = fPsi(4)*r*detJ
-
-    END SUBROUTINE areaQuad
-
-    !Computes element functions in point Xi
-    PURE FUNCTION fPsiQuad(self, Xi, nNodes) RESULT(fPsi)
-      IMPLICIT NONE
-
-      CLASS(meshCell2DCylQuad), INTENT(in):: self
-      REAL(8), INTENT(in):: Xi(1:3)
-      INTEGER, INTENT(in):: nNodes
-      REAL(8):: fPsi(1:nNodes)
-
-      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(self, Xi, nNodes) RESULT(dPsi)
-      IMPLICIT NONE
-
-      CLASS(meshCell2DCylQuad), INTENT(in):: self
-      REAL(8), INTENT(in):: Xi(1:3)
-      INTEGER, INTENT(in):: nNodes
-      REAL(8):: dPsi(1:3,1:nNodes)
-
-      dPsi = 0.D0
-
-      dPsi(1,:) = (/ -(1.D0 - Xi(2)), &
-                      (1.D0 - Xi(2)), &
-                      (1.D0 + Xi(2)), &
-                     -(1.D0 + Xi(2)) /)
-
-      dPsi(2,:) = (/ -(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 SUBROUTINE partialDerQuad(self, nNodes, dPsi, dz, dr)
       IMPLICIT NONE
 
       CLASS(meshCell2DCylQuad), INTENT(in):: self
       INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
+      INTEGER:: n(1:nNodes)
-      REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
 
-      dz = (/ DOT_PRODUCT(dPsi(1,1:4),self%z(1:4)), &
+      n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /)
-              DOT_PRODUCT(dPsi(2,1:4),self%z(1:4)) /)
-      dr = (/ DOT_PRODUCT(dPsi(1,1:4),self%r(1:4)), &
-              DOT_PRODUCT(dPsi(2,1:4),self%r(1:4)) /)
 
-    END SUBROUTINE partialDerQuad
+    END FUNCTION getNodesQuad
 
     !Random position in quadrilateral volume
     FUNCTION randPosCellQuad(self) RESULT(r)
@@ -410,74 +318,64 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION randPosCellQuad
 
-    !Computes element local stiffness matrix
+    !Computes element functions in point Xi
-    PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK)
+    PURE FUNCTION fPsiQuad(Xi, nNodes) RESULT(fPsi)
-      USE moduleConstParam, ONLY: PI2
+      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(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.D0 - Xi(2)), &
+                      (1.D0 - Xi(2)), &
+                      (1.D0 + Xi(2)), &
+                     -(1.D0 + Xi(2)) /)
+
+      dPsi(2,:) = (/ -(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
 
       CLASS(meshCell2DCylQuad), INTENT(in):: self
       INTEGER, INTENT(in):: nNodes
-      REAL(8):: localK(1:nNodes,1:nNodes)
+      REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
-      REAL(8):: Xi(1:3)
+      REAL(8):: pDer(1:3, 1:3)
-      REAL(8):: fPsi(1:4), dPsi(1:3,1:4)
-      REAL(8):: r
-      REAL(8):: invJ(1:3,1:3), detJ
-      INTEGER:: l, m
 
-      localK=0.D0
+      pDer = 0.D0
-      Xi=0.D0
-      !Start 2D Gauss Quad Integral
-      DO l=1, 3
-        Xi(2) = corQuad(l)
-        DO m = 1, 3
-          Xi(1)  = corQuad(m)
-          fPsi   = self%fPsi(Xi, 4)
-          dPsi   = self%dPsi(Xi, 4)
-          detJ   = self%detJac(Xi, 4, dPsi)
-          invJ   = self%invJac(Xi, 4, dPsi)
-          r      = DOT_PRODUCT(fPsi,self%r)
-          localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
-                                             MATMUL(invJ,dPsi))* &
-                            r*wQuad(l)*wQuad(m)/detJ
 
-        END DO
+      pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%z(1:4)), &
-      END DO
+                        DOT_PRODUCT(dPsi(2,1:4),self%z(1:4)) /)
-      localK = localK*PI2
+      pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%r(1:4)), &
+                        DOT_PRODUCT(dPsi(2,1:4),self%r(1:4)) /)
 
-    END FUNCTION elemKQuad
+    END FUNCTION partialDerQuad
-
-    !Computes the local source vector for a force f
-    PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF)
-      USE moduleConstParam, ONLY: PI2
-      IMPLICIT NONE
-
-      CLASS(meshCell2DCylQuad), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in):: source(1:nNodes)
-      REAL(8):: localF(1:nNodes)
-      REAL(8):: Xi(1:3)
-      REAL(8):: fPsi(1:4)
-      REAL(8):: r
-      REAL(8):: detJ, f
-      INTEGER:: l, m
-
-      localF = 0.D0
-      Xi     = 0.D0
-      DO l=1, 3
-        Xi(1) = corQuad(l)
-        DO m = 1, 3
-          Xi(2)  = corQuad(m)
-          detJ   = self%detJac(Xi, 4)
-          fPsi = self%fPsi(Xi, 4)
-          r      = DOT_PRODUCT(fPsi,self%r)
-          f      = DOT_PRODUCT(fPsi,source)
-          localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ
-
-         END DO
-      END DO
-      localF = localF*PI2
-
-    END FUNCTION elemFQuad
 
     PURE FUNCTION gatherEFQuad(self, Xi) RESULT(array)
       IMPLICIT NONE
@@ -521,6 +419,80 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION gatherMFQuad
 
+    !Computes element local stiffness matrix
+    PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK)
+      USE moduleConstParam, ONLY: PI2
+      IMPLICIT NONE
+
+      CLASS(meshCell2DCylQuad), INTENT(in):: self
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: localK(1:nNodes,1:nNodes)
+      REAL(8):: Xi(1:3)
+      REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
+      REAL(8):: pDer(1:3, 1:3)
+      REAL(8):: r
+      REAL(8):: invJ(1:3,1:3), detJ
+      INTEGER:: l, m
+
+      localK=0.D0
+      Xi=0.D0
+      !Start 2D Gauss Quad Integral
+      DO l=1, 3
+        Xi(2) = corQuad(l)
+        DO m = 1, 3
+          Xi(1)  = corQuad(m)
+          dPsi   = self%dPsi(Xi, 4)
+          pDer   = self%partialDer(4, dPsi)
+          detJ   = self%detJac(pDer)
+          invJ   = self%invJac(pDer)
+          fPsi   = self%fPsi(Xi, 4)
+          r      = DOT_PRODUCT(fPsi,self%r)
+          localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
+                                             MATMUL(invJ,dPsi))* &
+                            r*wQuad(l)*wQuad(m)/detJ
+
+        END DO
+      END DO
+      localK = localK*PI2
+
+    END FUNCTION elemKQuad
+
+    !Computes the local source vector for a force f
+    PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF)
+      USE moduleConstParam, ONLY: PI2
+      IMPLICIT NONE
+
+      CLASS(meshCell2DCylQuad), INTENT(in):: self
+      INTEGER, INTENT(in):: nNodes
+      REAL(8), INTENT(in):: source(1:nNodes)
+      REAL(8):: localF(1:nNodes)
+      REAL(8):: Xi(1:3)
+      REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
+      REAL(8):: pDer(1:3, 1:3)
+      REAL(8):: r
+      REAL(8):: detJ, f
+      INTEGER:: l, m
+
+      localF = 0.D0
+      Xi     = 0.D0
+      DO l=1, 3
+        Xi(1) = corQuad(l)
+        DO m = 1, 3
+          Xi(2)  = corQuad(m)
+          dPsi   = self%dPsi(Xi, 4)
+          pDer   = self%partialDer(4, dPsi)
+          detJ   = self%detJac(pDer)
+          fPsi = self%fPsi(Xi, 4)
+          r      = DOT_PRODUCT(fPsi,self%r)
+          f      = DOT_PRODUCT(fPsi,source)
+          localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ
+
+         END DO
+      END DO
+      localF = localF*PI2
+
+    END FUNCTION elemFQuad
+
     !Checks if a particle is inside a quad element
     PURE FUNCTION insideQuad(Xi) RESULT(ins)
       IMPLICIT NONE
@@ -533,18 +505,6 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION insideQuad
 
-    !Gets nodes from quadrilateral element
-    PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n)
-      IMPLICIT NONE
-
-      CLASS(meshCell2DCylQuad), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      INTEGER:: n(1:nNodes)
-
-      n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /)
-
-    END FUNCTION getNodesQuad
-
     !Transforms physical coordinates to element coordinates
     PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) 
       IMPLICIT NONE
@@ -554,6 +514,7 @@ MODULE moduleMesh2DCyl
       REAL(8):: Xi(1:3)
       REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
       REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
+      REAL(8):: pDer(1:3, 1:3)
       REAL(8):: conv
 
       !Iterative newton method to transform coordinates
@@ -562,8 +523,9 @@ MODULE moduleMesh2DCyl
 
       DO WHILE(conv > 1.D-2)
         dPsi = self%dPsi(XiO, 4)
-        invJ = self%invJac(XiO, 4, dPsi)
+        pDer = self%partialDer(4, dPsi)
-        detJ = self%detJac(XiO, 4, dPsi)
+        detJ = self%detJac(pDer)
+        invJ = self%invJac(pDer)
         fPsi = self%fPsi(XiO, 4)
         f = (/ DOT_PRODUCT(fPsi,self%z), &
                DOT_PRODUCT(fPsi,self%r), &
@@ -578,31 +540,69 @@ MODULE moduleMesh2DCyl
     END FUNCTION phy2logQuad
 
     !Gets the next element for a logical position Xi
-    SUBROUTINE nextElementQuad(self, Xi, nextElement)
+    SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement)
       IMPLICIT NONE
 
       CLASS(meshCell2DCylQuad), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
-      CLASS(meshElement), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
       REAL(8):: XiArray(1:4)
       INTEGER:: nextInt
 
       XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /)
       nextInt = MAXLOC(XiArray,1)
       !Selects the higher value of directions and searches in that direction
-      NULLIFY(nextElement)
+      NULLIFY(neighbourElement)
       SELECT CASE (nextInt)
       CASE (1)
-        nextElement => self%e1
+        neighbourElement => self%e1
       CASE (2)
-        nextElement => self%e2
+        neighbourElement => self%e2
       CASE (3)
-        nextElement => self%e3
+        neighbourElement => self%e3
       CASE (4)
-        nextElement => self%e4
+        neighbourElement => self%e4
       END SELECT
 
-    END SUBROUTINE nextElementQuad
+    END SUBROUTINE neighbourElementQuad
+
+    !Computes element area
+    PURE SUBROUTINE areaQuad(self)
+      USE moduleConstParam, ONLY: PI8
+      IMPLICIT NONE
+
+      CLASS(meshCell2DCylQuad), INTENT(inout):: self
+      REAL(8):: r, Xi(1:3)
+      REAL(8):: detJ
+      REAL(8):: fPsi(1:4)
+      REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3)
+
+      self%volume  = 0.D0
+      self%arNodes = 0.D0
+      !2D 1 point Gauss Quad Integral
+      Xi = 0.D0
+      dPsi = self%dPsi(Xi, 4)
+      pDer = self%partialDer(4, dPsi)
+      detJ = self%detJac(pDer)*PI8 !4*2*pi
+      fPsi = self%fPsi(Xi, 4)
+      !Computes total volume of the cell
+      r = DOT_PRODUCT(fPsi,self%r)
+      self%volume  = r*detJ
+      !Computes volume per node
+      Xi = (/-5.D-1, -5.D-1, 0.D0/)
+      r = self%gatherF(Xi, 4, self%r)
+      self%arNodes(1) = fPsi(1)*r*detJ
+      Xi = (/ 5.D-1, -5.D-1, 0.D0/)
+      r = self%gatherF(Xi, 4, self%r)
+      self%arNodes(2) = fPsi(2)*r*detJ
+      Xi = (/ 5.D-1,  5.D-1, 0.D0/)
+      r = self%gatherF(Xi, 4, self%r)
+      self%arNodes(3) = fPsi(3)*r*detJ
+      Xi = (/-5.D-1,  5.D-1, 0.D0/)
+      r = self%gatherF(Xi, 4, self%r)
+      self%arNodes(4) = fPsi(4)*r*detJ
+
+    END SUBROUTINE areaQuad
 
     !TRIA ELEMENT
     !Init tria element
@@ -645,6 +645,18 @@ MODULE moduleMesh2DCyl
 
     END SUBROUTINE initCellTria2DCyl
 
+    !Gets node indexes from triangular element
+    PURE FUNCTION getNodesTria(self, nNodes) RESULT(n)
+      IMPLICIT NONE
+
+      CLASS(meshCell2DCylTria), INTENT(in):: self
+      INTEGER, INTENT(in):: nNodes
+      INTEGER:: n(1:nNodes)
+
+      n = (/self%n1%n, self%n2%n, self%n3%n /)
+
+    END FUNCTION getNodesTria
+
     !Random position in quadrilateral volume
     FUNCTION randPosCellTria(self) RESULT(r)
       USE moduleRandom
@@ -667,36 +679,10 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION randPosCellTria
 
-    !Calculates area for triangular element
-    PURE SUBROUTINE areaTria(self)
-      USE moduleConstParam, ONLY: PI
-      IMPLICIT NONE
-
-      CLASS(meshCell2DCylTria), INTENT(inout):: self
-      REAL(8):: Xi(1:3)
-      REAL(8):: r
-      REAL(8):: detJ
-      REAL(8):: fPsi(1:3)
-
-      self%volume  = 0.D0
-      self%arNodes = 0.D0
-      !2D 1 point Gauss Quad Integral
-      Xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
-      detJ  = self%detJac(Xi, 3)*PI !2PI*1/2
-      fPsi = self%fPsi(Xi, 4)
-      !Computes total volume of the cell
-      r = DOT_PRODUCT(fPsi,self%r)
-      self%volume  =      r*detJ
-      !Computes volume per node
-      self%arNodes = fPsi*r*detJ
-
-    END SUBROUTINE areaTria
-
     !Shape functions for triangular element
-    PURE FUNCTION fPsiTria(self, Xi, nNodes) RESULT(fPsi)
+    PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi)
       IMPLICIT NONE
 
-      CLASS(meshCell2DCylTria), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
       INTEGER, INTENT(in):: nNodes
       REAL(8):: fPsi(1:nNodes)
@@ -708,10 +694,9 @@ MODULE moduleMesh2DCyl
     END FUNCTION fPsiTria
 
     !Derivative element function at coordinates Xi
-    PURE FUNCTION dPsiTria(self, Xi, nNodes) RESULT(dPsi)
+    PURE FUNCTION dPsiTria(Xi, nNodes) RESULT(dPsi)
       IMPLICIT NONE
 
-      CLASS(meshCell2DCylTria), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
       INTEGER, INTENT(in):: nNodes
       REAL(8):: dPsi(1:3,1:nNodes)
@@ -723,84 +708,22 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION dPsiTria
 
-    PURE SUBROUTINE partialDerTria(self, nNodes, dPsi, dz, dr)
+    PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer)
       IMPLICIT NONE
 
       CLASS(meshCell2DCylTria), INTENT(in):: self
       INTEGER, INTENT(in):: nNodes
       REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
-      REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
+      REAL(8):: pDer(1:3, 1:3)
 
-      dz = (/ DOT_PRODUCT(dPsi(1,:),self%z), &
+      pDer = 0.D0
-              DOT_PRODUCT(dPsi(2,:),self%z) /)
-      dr = (/ DOT_PRODUCT(dPsi(1,:),self%r), &
-              DOT_PRODUCT(dPsi(2,:),self%r) /)
 
-    END SUBROUTINE partialDerTria
+      pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%z(1:3)), &
+                        DOT_PRODUCT(dPsi(2,1:3),self%z(1:3)) /)
+      pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%r(1:3)), &
+                        DOT_PRODUCT(dPsi(2,1:3),self%r(1:3)) /)
 
-    !Computes element local stiffness matrix
+    END FUNCTION partialDerTria
-    PURE FUNCTION elemKTria(self, nNodes) RESULT(localK)
-      USE moduleConstParam, ONLY: PI2
-      IMPLICIT NONE
-
-      CLASS(meshCell2DCylTria), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      REAL(8):: localK(1:nNodes,1:nNodes)
-      REAL(8):: Xi(1:3)
-      REAL(8):: r
-      REAL(8):: fPsi(1:3), dPsi(1:3,1:3)
-      REAL(8):: invJ(1:3,1:3), detJ
-      INTEGER:: l
-
-      localK=0.D0
-      Xi=0.D0
-      !Start 2D Gauss Quad Integral
-      DO l=1, 4
-        Xi(1)  = Xi1Tria(l)
-        Xi(2)  = Xi2Tria(l)
-        dPsi   = self%dPsi(Xi, 4)
-        detJ   = self%detJac(Xi, 3, dPsi)
-        invJ   = self%invJac(Xi, 3, dPsi)
-        fPsi   = self%fPsi(Xi, 4)
-        r      = DOT_PRODUCT(fPsi,self%r)
-        localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ
-
-      END DO
-      localK = localK*PI2
-
-    END FUNCTION elemKTria
-
-    !Computes element local source vector
-    PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF)
-      USE moduleConstParam, ONLY: PI2
-      IMPLICIT NONE
-
-      CLASS(meshCell2DCylTria), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in):: source(1:nNodes)
-      REAL(8):: localF(1:nNodes)
-      REAL(8):: fPsi(1:3)
-      REAL(8):: Xi(1:3)
-      REAL(8):: r
-      REAL(8):: detJ, f
-      INTEGER:: l
-
-      localF = 0.D0
-      Xi     = 0.D0
-      !Start 2D Gauss Quad Integral
-      DO l=1, 4
-        Xi(1)  = Xi1Tria(l)
-        Xi(2)  = Xi2Tria(l)
-        detJ   = self%detJac(Xi, 3)
-        fPsi   = self%fPsi(Xi, 3)
-        r      = DOT_PRODUCT(fPsi,self%r)
-        f      = DOT_PRODUCT(fPsi,source)
-        localF = localF + r*f*fPsi*wTria(l)*detJ
-
-      END DO
-      localF = localF*PI2
-
-    END FUNCTION elemFTria
 
     PURE FUNCTION gatherEFTria(self, Xi) RESULT(array)
       IMPLICIT NONE
@@ -840,6 +763,75 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION gatherMFTria
 
+    !Computes element local stiffness matrix
+    PURE FUNCTION elemKTria(self, nNodes) RESULT(localK)
+      USE moduleConstParam, ONLY: PI2
+      IMPLICIT NONE
+
+      CLASS(meshCell2DCylTria), INTENT(in):: self
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: localK(1:nNodes,1:nNodes)
+      REAL(8):: Xi(1:3)
+      REAL(8):: r
+      REAL(8):: fPsi(1:3), dPsi(1:3,1:3)
+      REAL(8):: pDer(1:3, 1:3)
+      REAL(8):: invJ(1:3,1:3), detJ
+      INTEGER:: l
+
+      localK=0.D0
+      Xi=0.D0
+      !Start 2D Gauss Quad Integral
+      DO l=1, 4
+        Xi(1)  = Xi1Tria(l)
+        Xi(2)  = Xi2Tria(l)
+        dPsi   = self%dPsi(Xi, 3)
+        pDer   = self%partialDer(3, dPsi)
+        detJ   = self%detJac(pDer)
+        invJ   = self%invJac(pDer)
+        fPsi   = self%fPsi(Xi, 3)
+        r      = DOT_PRODUCT(fPsi,self%r)
+        localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ
+
+      END DO
+      localK = localK*PI2
+
+    END FUNCTION elemKTria
+
+    !Computes element local source vector
+    PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF)
+      USE moduleConstParam, ONLY: PI2
+      IMPLICIT NONE
+
+      CLASS(meshCell2DCylTria), INTENT(in):: self
+      INTEGER, INTENT(in):: nNodes
+      REAL(8), INTENT(in):: source(1:nNodes)
+      REAL(8):: localF(1:nNodes)
+      REAL(8):: fPsi(1:3)
+      REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3)
+      REAL(8):: Xi(1:3)
+      REAL(8):: r
+      REAL(8):: detJ, f
+      INTEGER:: l
+
+      localF = 0.D0
+      Xi     = 0.D0
+      !Start 2D Gauss Quad Integral
+      DO l=1, 4
+        Xi(1)  = Xi1Tria(l)
+        Xi(2)  = Xi2Tria(l)
+        dPsi   = self%dPsi(Xi, 3)
+        pDer   = self%partialDer(3, dPsi)
+        detJ   = self%detJac(pDer)
+        fPsi   = self%fPsi(Xi, 3)
+        r      = DOT_PRODUCT(fPsi,self%r)
+        f      = DOT_PRODUCT(fPsi,source)
+        localF = localF + r*f*fPsi*wTria(l)*detJ
+
+      END DO
+      localF = localF*PI2
+
+    END FUNCTION elemFTria
+
     PURE FUNCTION insideTria(Xi) RESULT(ins)
       IMPLICIT NONE
 
@@ -852,18 +844,6 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION insideTria
 
-    !Gets node indexes from triangular element
-    PURE FUNCTION getNodesTria(self, nNodes) RESULT(n)
-      IMPLICIT NONE
-
-      CLASS(meshCell2DCylTria), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      INTEGER:: n(1:nNodes)
-
-      n = (/self%n1%n, self%n2%n, self%n3%n /)
-
-    END FUNCTION getNodesTria
-
     !Transforms physical coordinates to element coordinates
     PURE FUNCTION phy2logTria(self,r) RESULT(Xi)
       IMPLICIT NONE
@@ -871,96 +851,97 @@ MODULE moduleMesh2DCyl
       CLASS(meshCell2DCylTria), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
-      REAL(8):: invJ(1:3,1:3), detJ
       REAL(8):: deltaR(1:3)
-      REAL(8):: dPsi(1:3,1:3)
+      REAL(8):: dPsi(1:3, 1:3)
+      REAL(8):: pDer(1:3, 1:3)
+      REAL(8):: invJ(1:3,  1:3), detJ
 
       !Direct method to convert coordinates
       Xi     = 0.D0
       deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /)
       dPsi   = self%dPsi(Xi, 3)
-      invJ   = self%invJac(Xi, 3,  dPsi)
+      pDer   = self%partialDer(3, dPsi)
-      detJ   = self%detJac(Xi, 3,  dPsi)
+      invJ   = self%invJac(pDer)
+      detJ   = self%detJac(pDer)
       Xi     = MATMUL(invJ,deltaR)/detJ
 
     END FUNCTION phy2logTria
 
-    SUBROUTINE nextElementTria(self, Xi, nextElement)
+    SUBROUTINE neighbourElementTria(self, Xi, neighbourElement)
       IMPLICIT NONE
 
       CLASS(meshCell2DCylTria), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
-      CLASS(meshElement), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
       REAL(8):: XiArray(1:3)
       INTEGER:: nextInt
 
       XiArray = (/ Xi(2), 1.D0-Xi(1)-Xi(2), Xi(1) /)
       nextInt = MINLOC(XiArray,1)
-      NULLIFY(nextElement)
+      NULLIFY(neighbourElement)
       SELECT CASE (nextInt)
       CASE (1)
-        nextElement => self%e1
+        neighbourElement => self%e1
       CASE (2)
-        nextElement => self%e2
+        neighbourElement => self%e2
       CASE (3)
-        nextElement => self%e3
+        neighbourElement => self%e3
       END SELECT
 
-    END SUBROUTINE nextElementTria
+    END SUBROUTINE neighbourElementTria
+
+    !Calculates area for triangular element
+    PURE SUBROUTINE areaTria(self)
+      USE moduleConstParam, ONLY: PI
+      IMPLICIT NONE
+
+      CLASS(meshCell2DCylTria), INTENT(inout):: self
+      REAL(8):: Xi(1:3)
+      REAL(8):: r
+      REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3)
+      REAL(8):: detJ
+      REAL(8):: fPsi(1:3)
+
+      self%volume  = 0.D0
+      self%arNodes = 0.D0
+      !2D 1 point Gauss Quad Integral
+      Xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
+      dPsi = self%dPsi(Xi, 3)
+      pDer = self%partialDer(3, dPsi)
+      detJ = self%detJac(pDer)*PI !2PI*1/2
+      fPsi = self%fPsi(Xi, 4)
+      !Computes total volume of the cell
+      r = DOT_PRODUCT(fPsi,self%r)
+      self%volume  =      r*detJ
+      !Computes volume per node
+      self%arNodes = fPsi*r*detJ
+
+    END SUBROUTINE areaTria
 
     !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS
     !Computes element Jacobian determinant
-    PURE FUNCTION detJ2DCyl(self, Xi, nNodes, dPsi_in) RESULT(dJ)
+    PURE FUNCTION detJ2DCyl(pDer) RESULT(dJ)
       IMPLICIT NONE
 
-      CLASS(meshCell2DCyl), INTENT(in):: self
+      REAL(8), INTENT(in):: pDer(1:3, 1:3)
-      REAL(8), INTENT(in):: Xi(1:3)
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
       REAL(8):: dJ
-      REAL(8):: dPsi(1:3,1:nNodes)
-      REAL(8):: dz(1:2), dr(1:2)
 
-      IF(PRESENT(dPsi_in)) THEN
+      dJ = pDer(1,1)*pDer(2,2)-pDer(1,2)*pDer(2,1)
-        dPsi = dPsi_in
-
-      ELSE
-        dPsi = self%dPsi(Xi, nNodes)
-
-      END IF
-
-      CALL self%partialDer(nNodes, dPsi, dz, dr)
-
-      dJ = dz(1)*dr(2)-dz(2)*dr(1)
 
     END FUNCTION detJ2DCyl
 
     !Computes element Jacobian inverse matrix (without determinant)
-    PURE FUNCTION invJ2DCyl(self, Xi, nNodes, dPsi_in) RESULT(invJ)
+    PURE FUNCTION invJ2DCyl(pDer) RESULT(invJ)
       IMPLICIT NONE
 
-      CLASS(meshCell2DCyl), INTENT(in):: self
+      REAL(8), INTENT(in):: pDer(1:3, 1:3)
-      REAL(8), INTENT(in):: Xi(1:3)
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
       REAL(8):: invJ(1:3,1:3)
-      REAL(8):: dPsi(1:3,1:nNodes)
-      REAL(8):: dz(1:2), dr(1:2)
-
-      IF(PRESENT(dPsi_in)) THEN
-        dPsi=dPsi_in
-
-      ELSE
-        dPsi = self%dPsi(Xi, 4)
-
-      END IF
 
       invJ = 0.D0
 
-      CALL self%partialDer(nNodes, dPsi, dz, dr)
+      invJ(1,1:2) = (/  pDer(2,2), -pDer(1,2) /)
-
+      invJ(2,1:2) = (/ -pDer(2,1),  pDer(1,1) /)
-      invJ(1,1:2) = (/  dr(2), -dz(2) /)
+      invJ(3,3)   = 1.D0
-      invJ(2,1:2) = (/ -dr(1),  dz(1) /)
 
     END FUNCTION invJ2DCyl
 

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -11,6 +11,7 @@ MODULE moduleMesh3DCart
     !Element coordinates
     REAL(8):: x, y, z
     CONTAINS
+      !meshNode DEFERRED PROCEDURES
       PROCEDURE, PASS:: init           => initNode3DCart
       PROCEDURE, PASS:: getCoordinates => getCoord3DCart
 
@@ -23,36 +24,18 @@ MODULE moduleMesh3DCart
     !Connectivity to nodes
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL()
     CONTAINS
+      !meshEdge DEFERRED PROCEDURES
       PROCEDURE, PASS::   init         => initEdge3DCartTria
       PROCEDURE, PASS::   getNodes     => getNodes3DCartTria
       PROCEDURE, PASS::   intersection => intersection3DCartTria
       PROCEDURE, PASS::   randPos      => randPosEdgeTria
+      !PARTICULAR PROCEDURES
       PROCEDURE, NOPASS:: fPsi         => fPsiEdgeTria
 
   END TYPE meshEdge3DCartTria
 
-  TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell3DCart
-    CONTAINS
-      PROCEDURE, PASS:: detJac => detJ3DCart
-      PROCEDURE, PASS:: invJac => invJ3DCart
-      PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
-  
-  END TYPE meshCell3DCart
-
-  ABSTRACT INTERFACE
-    PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx, dy, dz)
-      IMPORT meshCell3DCart
-      CLASS(meshCell3DCart), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
-      REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz
-
-    END SUBROUTINE partialDer_interface
-
-  END INTERFACE
-
   !Tetrahedron volume element
-  TYPE, PUBLIC, EXTENDS(meshCell3DCart):: meshCell3DCartTetra
+  TYPE, PUBLIC, EXTENDS(meshCell):: meshCell3DCartTetra
     !Element Coordinates
     REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0, z(1:4) = 0.D0
     !Connectivity to nodes
@@ -60,22 +43,24 @@ MODULE moduleMesh3DCart
     !Connectivity to adjacent elements
     CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
     CONTAINS
+      !meshCell DEFERRED PROCEDURES
       PROCEDURE, PASS::   init                => initCellTetra
+      PROCEDURE, PASS::   getNodes            => getNodesTetra
       PROCEDURE, PASS::   randPos             => randPosCellTetra
-      PROCEDURE, PASS::   calcCell     => volumeTetra
+      PROCEDURE, NOPASS:: fPsi                => fPsiTetra
-      PROCEDURE, PASS:: fPsi        => fPsiTetra
+      PROCEDURE, NOPASS:: dPsi                => dPsiTetra
-      PROCEDURE, PASS:: dPsi        => dPsiTetra
-      PROCEDURE, NOPASS, PRIVATE:: dPsiXi1 => dPsiTetraXi1
-      PROCEDURE, NOPASS, PRIVATE:: dPsiXi2 => dPsiTetraXi2
       PROCEDURE, PASS::   partialDer          => partialDerTetra
-      PROCEDURE, PASS::   elemK       => elemKTetra
+      PROCEDURE, NOPASS:: detJac              => detJ3DCart
-      PROCEDURE, PASS::   elemF       => elemFTetra
+      PROCEDURE, NOPASS:: invJac              => invJ3DCart
       PROCEDURE, PASS::   gatherElectricField => gatherEFTetra
       PROCEDURE, PASS::   gatherMagneticField => gatherMFTetra
+      PROCEDURE, PASS::   elemK               => elemKTetra
+      PROCEDURE, PASS::   elemF               => elemFTetra
       PROCEDURE, NOPASS:: inside              => insideTetra
-      PROCEDURE, PASS::   getNodes    => getNodesTetra
       PROCEDURE, PASS::   phy2log             => phy2logTetra
-      PROCEDURE, PASS::   nextElement => nextElementTetra
+      PROCEDURE, PASS::   neighbourElement    => neighbourElementTetra
+      !PARTICULAR PROCEDURES
+      PROCEDURE, PASS::   calcVol             => volumeTetra
       
   END TYPE meshCell3DCartTetra
 
@@ -227,9 +212,7 @@ MODULE moduleMesh3DCart
       IMPLICIT NONE
 
       REAL(8), INTENT(in):: Xi(1:3)
-      REAL(8), ALLOCATABLE:: fPsi(:)
+      REAL(8):: fPsi(1:3)
-
-      ALLOCATE(fPsi(1:3))
 
       fPsi(1) = 1.D0 - Xi(1) - Xi(2)
       fPsi(2) =        Xi(1)
@@ -268,7 +251,7 @@ MODULE moduleMesh3DCart
       self%z  = (/r1(3), r2(3), r3(3), r4(3)/)
 
       !Computes the element volume
-      CALL self%calcCell()
+      CALL self%calcVol()
 
       !Assign proportional volume to each node
       Xi = (/0.25D0, 0.25D0, 0.25D0/)
@@ -286,6 +269,17 @@ MODULE moduleMesh3DCart
 
     END SUBROUTINE initCellTetra
 
+    PURE FUNCTION getNodesTetra(self, nNodes) RESULT(n)
+      IMPLICIT NONE
+
+      CLASS(meshCell3DCartTetra), INTENT(in):: self
+      INTEGER, INTENT(in):: nNodes
+      INTEGER:: n(1:nNodes)
+
+      n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /)
+
+    END FUNCTION getNodesTetra
+
     !Random position in volume tetrahedron
     FUNCTION randPosCellTetra(self) RESULT(r)
       USE moduleRandom
@@ -308,24 +302,10 @@ MODULE moduleMesh3DCart
 
     END FUNCTION randPosCellTetra
 
-    !Computes the element volume
-    PURE SUBROUTINE volumeTetra(self)
-      IMPLICIT NONE
-
-      CLASS(meshCell3DCartTetra), INTENT(inout):: self
-      REAL(8):: Xi(1:3)
-
-      self%volume = 0.D0
-      Xi = (/0.25D0, 0.25D0, 0.25D0/)
-      self%volume = self%detJac(Xi, 4)
-
-    END SUBROUTINE volumeTetra
-
     !Computes element functions in point Xi
-    PURE FUNCTION fPsiTetra(self, Xi, nNodes) RESULT(fPsi)
+    PURE FUNCTION fPsiTetra(Xi, nNodes) RESULT(fPsi)
       IMPLICIT NONE
 
-      CLASS(meshCell3DCartTetra), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
       INTEGER, INTENT(in):: nNodes
       REAL(8):: fPsi(1:nNodes)
@@ -338,127 +318,45 @@ MODULE moduleMesh3DCart
     END FUNCTION fPsiTetra
 
     !Derivative element function at coordinates Xi
-    PURE FUNCTION dPsiTetra(self, Xi, nNodes) RESULT(dPsi)
+    PURE FUNCTION dPsiTetra(Xi, nNodes) RESULT(dPsi)
       IMPLICIT NONE
 
-      CLASS(meshCell3DCartTetra), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
       INTEGER, INTENT(in):: nNodes
       REAL(8):: dPsi(1:3, 1:nNodes)
 
       dPsi = 0.D0
 
-      dPsi(1,:) = dPsiTetraXi1(Xi(2), Xi(3))
+      dPsi(1,1:4) = (/ -1.D0, 1.D0, 0.D0, 0.D0 /)
-      dPsi(2,:) = dPsiTetraXi2(Xi(1), Xi(3))
+      dPsi(2,1:4) = (/ -1.D0, 0.D0, 1.D0, 0.D0 /)
-      dPsi(3,:) = dPsiTetraXi3(Xi(1), Xi(2))
+      dPsi(3,1:4) = (/ -1.D0, 0.D0, 0.D0, 1.D0 /)
 
     END FUNCTION dPsiTetra
 
-    !Derivative element function respect to Xi1
-    PURE FUNCTION dPsiTetraXi1(Xi2, Xi3) RESULT(dPsiXi1)
-      IMPLICIT NONE
-      REAL(8), INTENT(in):: Xi2, Xi3
-      REAL(8):: dPsiXi1(1:4)
-
-      dPsiXi1(1) = -1.D0
-      dPsiXi1(2) =  1.D0
-      dPsiXi1(3) =  0.D0
-      dPsiXi1(4) =  0.D0
-
-    END FUNCTION dPsiTetraXi1
-      
-    !Derivative element function respect to Xi2
-    PURE FUNCTION dPsiTetraXi2(Xi1, Xi3) RESULT(dPsiXi2)
-      IMPLICIT NONE
-      REAL(8), INTENT(in):: Xi1, Xi3
-      REAL(8):: dPsiXi2(1:4)
-
-      dPsiXi2(1) = -1.D0
-      dPsiXi2(2) =  0.D0
-      dPsiXi2(3) =  1.D0
-      dPsiXi2(4) =  0.D0
-
-    END FUNCTION dPsiTetraXi2
-      
-    !Derivative element function respect to Xi3
-    PURE FUNCTION dPsiTetraXi3(Xi1, Xi2) RESULT(dPsiXi3)
-      IMPLICIT NONE
-      REAL(8), INTENT(in):: Xi1, Xi2
-      REAL(8):: dPsiXi3(1:4)
-
-      dPsiXi3(1) = -1.D0
-      dPsiXi3(2) =  0.D0
-      dPsiXi3(3) =  0.D0
-      dPsiXi3(4) =  1.D0
-
-    END FUNCTION dPsiTetraXi3
-    
     !Computes the derivatives in global coordinates
-    PURE SUBROUTINE partialDerTetra(self, nNodes, dPsi, dx, dy, dz)
+    PURE FUNCTION partialDerTetra(self, nNodes, dPsi) RESULT(pDer)
       IMPLICIT NONE
 
       CLASS(meshCell3DCartTetra), INTENT(in):: self
       INTEGER, INTENT(in):: nNodes
       REAL(8), INTENT(in):: dPsi(1:3, 1:nNodes)
-      REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz
+      REAL(8):: pDer(1:3, 1:3)
 
-      dx(1) = DOT_PRODUCT(dPsi(1,:), self%x)
+      pDer = 0.D0
-      dx(2) = DOT_PRODUCT(dPsi(2,:), self%x)
-      dx(3) = DOT_PRODUCT(dPsi(3,:), self%x)
 
-      dy(1) = DOT_PRODUCT(dPsi(1,:), self%y)
+      pDer(1, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%x(1:4)), &
-      dy(2) = DOT_PRODUCT(dPsi(2,:), self%y)
+                        DOT_PRODUCT(dPsi(2,1:4), self%x(1:4)), &
-      dy(3) = DOT_PRODUCT(dPsi(3,:), self%y)
+                        DOT_PRODUCT(dPsi(3,1:4), self%x(1:4)) /)
 
-      dz(1) = DOT_PRODUCT(dPsi(1,:), self%z)
+      pDer(2, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%y(1:4)), &
-      dz(2) = DOT_PRODUCT(dPsi(2,:), self%z)
+                        DOT_PRODUCT(dPsi(2,1:4), self%y(1:4)), &
-      dz(3) = DOT_PRODUCT(dPsi(3,:), self%z)
+                        DOT_PRODUCT(dPsi(3,1:4), self%y(1:4)) /)
 
-    END SUBROUTINE partialDerTetra
+      pDer(3, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%z(1:4)), &
+                        DOT_PRODUCT(dPsi(2,1:4), self%z(1:4)), &
+                        DOT_PRODUCT(dPsi(3,1:4), self%z(1:4)) /)
 
-    PURE FUNCTION elemKTetra(self, nNodes) RESULT(localK)
+    END FUNCTION partialDerTetra
-      IMPLICIT NONE
-
-      CLASS(meshCell3DCartTetra), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      REAL(8):: localK(1:nNodes,1:nNodes)
-      REAL(8):: Xi(1:3)
-      REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
-      REAL(8):: invJ(1:3,1:3), detJ
-
-      localK = 0.D0
-      Xi    = 0.D0
-      !TODO: One point Gauss integral. Upgrade when possible
-      Xi    = (/ 0.25D0, 0.25D0, 0.25D0 /)
-      dPsi   = self%dPsi(Xi, 4)
-      detJ   = self%detJac(Xi, 4, dPsi)
-      invJ   = self%invJac(Xi, 4, dPsi)
-      fPsi = self%fPsi(Xi, 4)
-      localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ
-
-    END FUNCTION elemKTetra
-
-    PURE FUNCTION elemFTetra(self, nNodes, source) RESULT(localF)
-      IMPLICIT NONE
-
-      CLASS(meshCell3DCartTetra), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in):: source(1:nNodes)
-      REAL(8):: localF(1:nNodes)
-      REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
-      REAL(8):: Xi(1:3)
-      REAL(8):: detJ, f
-
-      localF = 0.D0
-      Xi   = 0.D0
-      Xi   = (/ 0.25D0, 0.25D0, 0.25D0 /)
-      dPsi = self%dPsi(Xi, 4)
-      detJ = self%detJac(Xi, 4, dPsi)
-      fPsi = self%fPsi(Xi, 4)
-      f    = DOT_PRODUCT(fPsi, source)
-      localF = f*fPsi*1.D0*detJ
-
-    END FUNCTION elemFTetra
 
     PURE FUNCTION gatherEFTetra(self, Xi) RESULT(array)
       IMPLICIT NONE
@@ -502,6 +400,54 @@ MODULE moduleMesh3DCart
 
     END FUNCTION gatherMFTetra
 
+    PURE FUNCTION elemKTetra(self, nNodes) RESULT(localK)
+      IMPLICIT NONE
+
+      CLASS(meshCell3DCartTetra), INTENT(in):: self
+      INTEGER, INTENT(in):: nNodes
+      REAL(8):: localK(1:nNodes,1:nNodes)
+      REAL(8):: Xi(1:3)
+      REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
+      REAL(8):: pDer(1:3, 1:3)
+      REAL(8):: invJ(1:3,1:3), detJ
+
+      localK = 0.D0
+      Xi    = 0.D0
+      !TODO: One point Gauss integral. Upgrade when possible
+      Xi     = (/ 0.25D0, 0.25D0, 0.25D0 /)
+      dPsi   = self%dPsi(Xi, 4)
+      pDer   = self%partialDer(4, dPsi)
+      detJ   = self%detJac(pDer)
+      invJ   = self%invJac(pDer)
+      fPsi = self%fPsi(Xi, 4)
+      localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ
+
+    END FUNCTION elemKTetra
+
+    PURE FUNCTION elemFTetra(self, nNodes, source) RESULT(localF)
+      IMPLICIT NONE
+
+      CLASS(meshCell3DCartTetra), INTENT(in):: self
+      INTEGER, INTENT(in):: nNodes
+      REAL(8), INTENT(in):: source(1:nNodes)
+      REAL(8):: localF(1:nNodes)
+      REAL(8):: Xi(1:3)
+      REAL(8):: fPsi(1:4), dPsi(1:3, 1:4)
+      REAL(8):: pDer(1:3, 1:3)
+      REAL(8):: detJ, f
+
+      localF = 0.D0
+      Xi   = 0.D0
+      Xi   = (/ 0.25D0, 0.25D0, 0.25D0 /)
+      dPsi = self%dPsi(Xi, 4)
+      pDer = self%partialDer(4, dPsi)
+      detJ = self%detJac(pDer)
+      fPsi = self%fPsi(Xi, 4)
+      f    = DOT_PRODUCT(fPsi, source)
+      localF = f*fPsi*1.D0*detJ
+
+    END FUNCTION elemFTetra
+
     PURE FUNCTION insideTetra(Xi) RESULT(ins)
       IMPLICIT NONE
 
@@ -515,121 +461,101 @@ MODULE moduleMesh3DCart
       
     END FUNCTION insideTetra
 
-    PURE FUNCTION getNodesTetra(self, nNodes) RESULT(n)
-      IMPLICIT NONE
-
-      CLASS(meshCell3DCartTetra), INTENT(in):: self
-      INTEGER, INTENT(in):: nNodes
-      INTEGER:: n(1:nNodes)
-
-      n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /)
-
-    END FUNCTION getNodesTetra
-
     PURE FUNCTION phy2logTetra(self,r) RESULT(Xi)
       IMPLICIT NONE
 
       CLASS(meshCell3DCartTetra), INTENT(in):: self
       REAL(8), INTENT(in):: r(1:3)
       REAL(8):: Xi(1:3)
+      REAL(8):: dPsi(1:3, 1:4)
+      REAL(8):: pDer(1:3, 1:3)
       REAL(8):: invJ(1:3, 1:3), detJ
       REAL(8):: deltaR(1:3)
-      REAL(8):: dPsi(1:3, 1:4)
 
       Xi     = 0.D0
       deltaR = (/r(1) - self%x(1), r(2) - self%y(1), r(3) - self%z(1) /)
       dPsi   = self%dPsi(Xi, 4)
-      invJ   = self%invJac(Xi, 4, dPsi)
+      pDer   = self%partialDer(4, dPsi)
-      detJ   = self%detJac(Xi, 4, dPsi)
+      invJ   = self%invJac(pDer)
+      detJ   = self%detJac(pDer)
       Xi     = MATMUL(invJ, deltaR)/detJ
 
     END FUNCTION phy2logTetra
     
-    SUBROUTINE nextElementTetra(self, Xi, nextElement)
+    SUBROUTINE neighbourElementTetra(self, Xi, neighbourElement)
       IMPLICIT NONE
 
       CLASS(meshCell3DCartTetra), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
-      CLASS(meshElement), POINTER, INTENT(out):: nextElement
+      CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
       REAL(8):: XiArray(1:4)
       INTEGER:: nextInt
 
       !TODO: Review when connectivity
       XiArray = (/ Xi(3), 1.D0 - Xi(1) - Xi(2) - Xi(3), Xi(2), Xi(1) /)
       nextInt = MINLOC(XiArray, 1)
-      NULLIFY(nextElement)
+      NULLIFY(neighbourElement)
       SELECT CASE(nextInt)
       CASE (1)
-        nextElement => self%e1
+        neighbourElement => self%e1
       CASE (2)
-        nextElement => self%e2
+        neighbourElement => self%e2
       CASE (3)
-        nextElement => self%e3
+        neighbourElement => self%e3
       CASE (4)
-        nextElement => self%e4
+        neighbourElement => self%e4
       END SELECT
 
-    END SUBROUTINE nextElementTetra
+    END SUBROUTINE neighbourElementTetra
+
+    !Computes the element volume
+    PURE SUBROUTINE volumeTetra(self)
+      IMPLICIT NONE
+
+      CLASS(meshCell3DCartTetra), INTENT(inout):: self
+      REAL(8):: Xi(1:3)
+      REAL(8):: dPsi(1:3, 1:4)
+      REAL(8):: pDer(1:3, 1:3)
+
+      self%volume = 0.D0
+      Xi = (/0.25D0, 0.25D0, 0.25D0/)
+      dPsi = self%dPsi(Xi, 4)
+      pDer = self%partialDer(4, dPsi)
+      self%volume = self%detJac(pDer)
+
+    END SUBROUTINE volumeTetra
 
     !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 3D
     !Computes element Jacobian determinant
-    PURE FUNCTION detJ3DCart(self, Xi, nNodes, dPsi_in) RESULT(dJ)
+    PURE FUNCTION detJ3DCart(pDer) RESULT(dJ)
       IMPLICIT NONE
 
-      CLASS(meshCell3DCart), INTENT(in)::self
+      REAL(8), INTENT(in):: pDer(1:3, 1:3)
-      REAL(8), INTENT(in):: Xi(1:3)
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3, 1:nNodes)
       REAL(8):: dJ
-      REAL(8):: dPsi(1:3, 1:nNodes)
-      REAL(8):: dx(1:3), dy(1:3), dz(1:3)
 
-      IF (PRESENT(dPsi_in)) THEN
+      dJ =   pDer(1,1)*(pDer(2,2)*pDer(3,3) - pDer(2,3)*pDer(3,2)) &
-        dPsi = dPsi_in
+           - pDer(1,2)*(pDer(2,1)*pDer(3,3) - pDer(2,3)*pDer(3,1)) &
-
+           + pDer(1,3)*(pDer(2,1)*pDer(3,2) - pDer(2,2)*pDer(3,1))
-      ELSE
-        dPsi = self%dPsi(Xi, 4)
-
-      END IF
-
-      CALL self%partialDer(nNodes, dPsi, dx, dy, dz)
-      dJ =   dx(1)*(dy(2)*dz(3) - dy(3)*dz(2)) &
-           - dx(2)*(dy(1)*dz(3) - dy(3)*dz(1)) &
-           + dx(3)*(dy(1)*dz(2) - dy(2)*dz(1))
 
     END FUNCTION detJ3DCart
 
-    PURE FUNCTION invJ3DCart(self, Xi, nNodes, dPsi_in) RESULT(invJ)
+    PURE FUNCTION invJ3DCart(pDer) RESULT(invJ)
       IMPLICIT NONE
 
-      CLASS(meshCell3DCart), INTENT(in):: self
+      REAL(8), INTENT(in):: pDer(1:3, 1:3)
-      REAL(8), INTENT(in):: Xi(1:3)
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3, 1:nNodes)
-      REAL(8):: dPsi(1:3, 1:nNodes)
-      REAL(8), DIMENSION(1:3):: dx, dy, dz
       REAL(8):: invJ(1:3,1:3)
 
-      IF(PRESENT(dPsi_in)) THEN
+      invJ(1,1:3) =  (/ (pDer(2,2)*pDer(3,3) - pDer(2,3)*pDer(3,2)), &
-        dPsi=dPsi_in
+                       -(pDer(2,1)*pDer(3,3) - pDer(2,3)*pDer(3,1)), &
+                        (pDer(2,1)*pDer(3,2) - pDer(2,2)*pDer(3,1)) /)
 
-      ELSE
+      invJ(2,1:3) = (/ -(pDer(1,2)*pDer(3,3) - pDer(1,3)*pDer(3,2)), &
-        dPsi = self%dPsi(Xi, 4)
+                        (pDer(1,1)*pDer(3,3) - pDer(1,3)*pDer(3,1)), &
+                       -(pDer(1,1)*pDer(3,2) - pDer(1,2)*pDer(3,1)) /)
 
-      END IF
+      invJ(3,1:3) =  (/ (pDer(1,2)*pDer(2,3) - pDer(1,3)*pDer(2,2)), &
-
+                       -(pDer(1,1)*pDer(2,3) - pDer(1,3)*pDer(2,1)), &
-      CALL self%partialDer(nNodes, dPsi, dx, dy, dz)
+                        (pDer(1,1)*pDer(2,2) - pDer(1,2)*pDer(2,1)) /)
-      invJ(1,1) =  (dy(2)*dz(3) - dy(3)*dz(2))
-      invJ(1,2) = -(dy(1)*dz(3) - dy(3)*dz(1))
-      invJ(1,3) =  (dy(1)*dz(2) - dy(2)*dz(1))
-
-      invJ(2,1) = -(dx(2)*dz(3) - dx(3)*dz(2))
-      invJ(2,2) =  (dx(1)*dz(3) - dx(3)*dz(1))
-      invJ(2,3) = -(dx(1)*dz(2) - dx(2)*dz(1))
-
-      invJ(3,1) =  (dx(2)*dy(3) - dx(3)*dy(2))
-      invJ(3,2) = -(dx(1)*dy(3) - dx(3)*dy(1))
-      invJ(3,3) =  (dx(1)*dy(2) - dx(2)*dy(1))
 
       invJ = TRANSPOSE(invJ)
 

src/modules/mesh/moduleMesh.f90

@@ -24,8 +24,10 @@ MODULE moduleMesh
     !Lock indicator for scattering
     INTEGER(KIND=OMP_LOCK_KIND):: lock
     CONTAINS
+      !DEFERED PROCEDURES
       PROCEDURE(initNode_interface),    DEFERRED, PASS:: init
       PROCEDURE(getCoord_interface),    DEFERRED, PASS:: getCoordinates
+      !GENERIC PROCEDURES
       PROCEDURE,                                  PASS:: resetOutput
 
   END TYPE meshNode
@@ -83,6 +85,7 @@ MODULE moduleMesh
     !Physical surface for the edge
     INTEGER:: physicalSurface
     CONTAINS
+      !DEFERED PROCEDURES
       PROCEDURE(initEdge_interface),         DEFERRED, PASS:: init
       PROCEDURE(getNodesEdge_interface),     DEFERRED, PASS:: getNodes
       PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection
@@ -166,37 +169,41 @@ MODULE moduleMesh
     !Total weight of particles inside cell
     REAL(8), ALLOCATABLE:: totalWeight(:)
     CONTAINS
+      !DEFERRED PROCEDURES
       !Init the cell
       PROCEDURE(initCell_interface),         DEFERRED, PASS::   init
       !Get the index of the nodes
       PROCEDURE(getNodesCell_interface),     DEFERRED, PASS::   getNodes
       !Calculate random position on the cell
-      PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos
+      PROCEDURE(randPosCell_interface),      DEFERRED, PASS::   randPos
       !Obtain functions and values of cell natural functions
-      PROCEDURE(fPsi_interface),   DEFERRED, PASS:: fPsi
+      PROCEDURE(fPsi_interface),             DEFERRED, NOPASS:: fPsi
-      PROCEDURE(dPsi_interface),   DEFERRED, PASS:: dPsi
+      PROCEDURE(dPsi_interface),             DEFERRED, NOPASS:: dPsi
-      PROCEDURE(detJac_interface), DEFERRED, PASS:: detJac
+      PROCEDURE(partialDer_interface),       DEFERRED, PASS::   partialDer
-      PROCEDURE(invJac_interface), DEFERRED, PASS:: invJac
+      PROCEDURE(detJac_interface),           DEFERRED, NOPASS:: detJac
-      !Scatter properties of particles on cell nodes
+      PROCEDURE(invJac_interface),           DEFERRED, NOPASS:: invJac
-      PROCEDURE, PASS:: scatter
-      !Gather value and spatial derivative on the nodes at position Xi
-      PROCEDURE, PASS, PRIVATE:: gatherF_scalar
-      PROCEDURE, PASS, PRIVATE:: gatherF_array
-      PROCEDURE, PASS, PRIVATE:: gatherDF_scalar
-      GENERIC:: gatherF  => gatherF_scalar, gatherF_array
-      GENERIC:: gatherDF => gatherDF_scalar
       !Procedures to get specific values in the node
       PROCEDURE(gatherArray_interface),      DEFERRED, PASS::   gatherElectricField
       PROCEDURE(gatherArray_interface),      DEFERRED, PASS::   gatherMagneticField
       !Compute K and F to solve PDE on the mesh
       PROCEDURE(elemK_interface),            DEFERRED, PASS::   elemK
       PROCEDURE(elemF_interface),            DEFERRED, PASS::   elemF
-      !Subroutines to find in which cell a particle is located
+      !Check if particle is inside the cell
-      PROCEDURE,                                  PASS::   findCell
       PROCEDURE(inside_interface),           DEFERRED, NOPASS:: inside
-      PROCEDURE(nextElement_interface), DEFERRED, PASS::   nextElement
       !Convert physical coordinates (r) into logical coordinates (Xi)
       PROCEDURE(phy2log_interface),          DEFERRED, PASS::   phy2log
+      !Returns the neighbour element based on particle position outside the cell
+      PROCEDURE(neighbourElement_interface), DEFERRED, PASS::   neighbourElement
+      !Scatter properties of particles on cell nodes
+      PROCEDURE, PASS:: scatter
+      !Subroutine to find in which cell a particle is located
+      PROCEDURE, PASS:: findCell
+      !Gather value and spatial derivative on the nodes at position Xi
+      PROCEDURE, PASS, PRIVATE:: gatherF_scalar
+      PROCEDURE, PASS, PRIVATE:: gatherF_array
+      PROCEDURE, PASS, PRIVATE:: gatherDF_scalar
+      GENERIC:: gatherF  => gatherF_scalar, gatherF_array
+      GENERIC:: gatherDF => gatherDF_scalar
 
   END TYPE meshCell
 
@@ -219,40 +226,44 @@ MODULE moduleMesh
 
     END FUNCTION getNodesCell_interface
 
-    PURE FUNCTION fPsi_interface(self, Xi, nNodes) RESULT(fPsi)
+    FUNCTION randPosCell_interface(self) RESULT(r)
       IMPORT:: meshCell
       CLASS(meshCell), INTENT(in):: self
+      REAL(8):: r(1:3)
+
+    END FUNCTION randPosCell_interface
+
+    PURE FUNCTION fPsi_interface(Xi, nNodes) RESULT(fPsi)
       REAL(8), INTENT(in):: Xi(1:3)
       INTEGER, INTENT(in):: nNodes
       REAL(8):: fPsi(1:nNodes)
 
     END FUNCTION fPsi_interface
 
-    PURE FUNCTION dPsi_interface(self, Xi, nNodes) RESULT(dPsi)
+    PURE FUNCTION dPsi_interface(Xi, nNodes) RESULT(dPsi)
-      IMPORT:: meshCell
-      CLASS(meshCell), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
       INTEGER, INTENT(in):: nNodes
       REAL(8):: dPsi(1:3, 1:nNodes)
 
     END FUNCTION dPsi_interface
 
-    PURE FUNCTION detJac_interface(self, Xi, nNodes, dPsi_in) RESULT(dJ)
+    PURE FUNCTION partialDer_interface(self, nNodes, dPsi) RESULT(pDer)
       IMPORT:: meshCell
       CLASS(meshCell), INTENT(in):: self
-      REAL(8), INTENT(in):: Xi(1:3)
       INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
+      REAL(8), INTENT(in):: dPsi(1:3, 1:nNodes)
+      REAL(8):: pDer(1:3, 1:3)
+
+    END FUNCTION partialDer_interface
+
+    PURE FUNCTION detJac_interface(pDer) RESULT(dJ)
+      REAL(8), INTENT(in):: pDer(1:3,1:3)
       REAL(8):: dJ
 
     END FUNCTION detJac_interface
 
-    PURE FUNCTION invJac_interface(self, Xi, nNodes, dPsi_in) RESULT(invJ)
+    PURE FUNCTION invJac_interface(pDer) RESULT(invJ)
-      IMPORT:: meshCell
+      REAL(8), INTENT(in):: pDer(1:3,1:3)
-      CLASS(meshCell), INTENT(in):: self
-      REAL(8), INTENT(in):: Xi(1:3)
-      INTEGER, INTENT(in):: nNodes
-      REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes)
       REAL(8):: invJ(1:3,1:3)
 
     END FUNCTION invJac_interface
@@ -282,13 +293,12 @@ MODULE moduleMesh
 
     END FUNCTION elemF_interface
 
-    SUBROUTINE nextElement_interface(self, Xi, nextElement)
+    PURE FUNCTION inside_interface(Xi) RESULT(ins)
-      IMPORT:: meshCell, meshElement
+      IMPORT:: meshCell
-      CLASS(meshCell), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
-      CLASS(meshElement), POINTER, INTENT(out):: nextElement
+      LOGICAL:: ins
 
-    END SUBROUTINE nextElement_interface
+    END FUNCTION inside_interface
 
     PURE FUNCTION phy2log_interface(self,r) RESULT(Xi) 
       IMPORT:: meshCell
@@ -298,19 +308,13 @@ MODULE moduleMesh
 
     END FUNCTION phy2log_interface
 
-    PURE FUNCTION inside_interface(Xi) RESULT(ins)
+    SUBROUTINE neighbourElement_interface(self, Xi, neighbourElement)
-      IMPORT:: meshCell
+      IMPORT:: meshCell, meshElement
-      REAL(8), INTENT(in):: Xi(1:3)
-      LOGICAL:: ins
-
-    END FUNCTION inside_interface
-
-    FUNCTION randPosVol_interface(self) RESULT(r)
-      IMPORT:: meshCell
       CLASS(meshCell), INTENT(in):: self
-      REAL(8):: r(1:3)
+      REAL(8), INTENT(in):: Xi(1:3)
+      CLASS(meshElement), POINTER, INTENT(out):: neighbourElement
 
-    END FUNCTION randPosVol_interface
+    END SUBROUTINE neighbourElement_interface
 
   END INTERFACE
 
@@ -332,11 +336,13 @@ MODULE moduleMesh
     TYPE(meshNodeCont), ALLOCATABLE:: nodes(:)
     !Array of cell elements
     TYPE(meshCellCont),  ALLOCATABLE:: cells(:)
+    !PROCEDURES SPECIFIC OF FILE TYPE
     PROCEDURE(readMesh_interface),     POINTER, PASS::   readMesh     => NULL()
     PROCEDURE(readInitial_interface),  POINTER, NOPASS:: readInitial  => NULL()
     PROCEDURE(connectMesh_interface),  POINTER, PASS::   connectMesh  => NULL()
     PROCEDURE(printColl_interface),    POINTER, PASS::   printColl    => NULL()
     CONTAINS
+      !GENERIC PROCEDURES
       PROCEDURE, PASS:: doCollisions
 
   END TYPE meshGeneric
@@ -345,7 +351,6 @@ MODULE moduleMesh
     !Reads the mesh from a file
     SUBROUTINE readMesh_interface(self, filename)
       IMPORT meshGeneric
-
       CLASS(meshGeneric), INTENT(inout):: self
       CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
 
@@ -363,7 +368,6 @@ MODULE moduleMesh
     !Connects cell and edges to the mesh
     SUBROUTINE connectMesh_interface(self)
       IMPORT meshGeneric
-
       CLASS(meshGeneric), INTENT(inout):: self
 
     END SUBROUTINE connectMesh_interface
@@ -371,7 +375,6 @@ MODULE moduleMesh
     !Prints number of collisions in each cell
     SUBROUTINE printColl_interface(self, t)
       IMPORT meshGeneric
-
       CLASS(meshGeneric), INTENT(inout):: self
       INTEGER, INTENT(in):: t
 
@@ -388,28 +391,21 @@ MODULE moduleMesh
     REAL(8), ALLOCATABLE, DIMENSION(:,:):: K 
     !Permutation matrix for P L U factorization
     INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV
+    !PROCEDURES SPECIFIC OF FILE TYPE
     PROCEDURE(printOutput_interface),  POINTER, PASS:: printOutput  => NULL()
     PROCEDURE(printEM_interface),      POINTER, PASS:: printEM      => NULL()
     PROCEDURE(doCoulomb_interface),    POINTER, PASS:: doCoulomb    => NULL()
     PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL()
     CONTAINS
+      !GENERIC PROCEDURES
       PROCEDURE, PASS:: constructGlobalK
 
   END TYPE meshParticles
 
   ABSTRACT INTERFACE
-    !Perform Coulomb Scattering
-    SUBROUTINE doCoulomb_interface(self)
-      IMPORT meshParticles
-
-      CLASS(meshParticles), INTENT(inout):: self
-
-    END SUBROUTINE doCoulomb_interface
-    
     !Prints Species data
     SUBROUTINE printOutput_interface(self, t)
       IMPORT meshParticles
-
       CLASS(meshParticles), INTENT(in):: self
       INTEGER, INTENT(in):: t
 
@@ -418,21 +414,25 @@ MODULE moduleMesh
     !Prints EM info
     SUBROUTINE printEM_interface(self, t)
       IMPORT meshParticles
-
       CLASS(meshParticles), INTENT(in):: self
       INTEGER, INTENT(in):: t
 
     END SUBROUTINE printEM_interface
 
+    !Perform Coulomb Scattering
+    SUBROUTINE doCoulomb_interface(self)
+      IMPORT meshParticles
+      CLASS(meshParticles), INTENT(inout):: self
+
+    END SUBROUTINE doCoulomb_interface
+    
     !Prints average values
     SUBROUTINE printAverage_interface(self)
       IMPORT meshParticles
-
       CLASS(meshParticles), INTENT(in):: self
 
     END SUBROUTINE printAverage_interface
 
-
   END INTERFACE
 
   TYPE(meshParticles), TARGET:: mesh
@@ -440,6 +440,7 @@ MODULE moduleMesh
   !Collision (MCC) mesh
   TYPE, EXTENDS(meshGeneric):: meshCollisions
     CONTAINS
+      !GENERIC PROCEDURES
 
   END TYPE meshCollisions
 
@@ -448,7 +449,6 @@ MODULE moduleMesh
   ABSTRACT INTERFACE
     SUBROUTINE readMeshColl_interface(self, filename)
       IMPORT meshCollisions
-
       CLASS(meshCollisions), INTENT(inout):: self
       CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
 
@@ -577,12 +577,14 @@ MODULE moduleMesh
       REAL(8), INTENT(in):: valNodes(1:nNodes)
       REAL(8):: df(1:3)
       REAL(8):: dPsi(1:3, 1:nNodes)
+      REAL(8):: pDer(1:3,1:3)
       REAL(8):: dPsiR(1:3, 1:nNodes)
       REAL(8):: invJ(1:3, 1:3), detJ
 
       dPsi  = self%dPsi(Xi, nNodes)
-      detJ  = self%detJac(Xi, nNodes, dPsi)
+      pDer  = self%partialDer(nNodes, dPsi)
-      invJ  = self%invJac(Xi, nNodes, dPsi)
+      detJ  = self%detJac(pDer)
+      invJ  = self%invJac(pDer)
       dPsiR =  MATMUL(invJ, dPsi)/detJ
       df    = (/ DOT_PRODUCT(dPsiR(1,:), valNodes), &
                  DOT_PRODUCT(dPsiR(2,:), valNodes), &
@@ -637,7 +639,7 @@ MODULE moduleMesh
       CLASS(particle), INTENT(inout), TARGET:: part
       CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell
       REAL(8):: Xi(1:3)
-      CLASS(meshElement), POINTER:: nextElement
+      CLASS(meshElement), POINTER:: neighbourElement
       INTEGER:: sp
 
       Xi = self%phy2log(part%r)
@@ -655,16 +657,16 @@ MODULE moduleMesh
 
       ELSE
         !If not, searches for a neighbour and repeats the process.
-        CALL self%nextElement(Xi, nextElement)
+        CALL self%neighbourElement(Xi, neighbourElement)
         !Defines the next step
-        SELECT TYPE(nextElement)
+        SELECT TYPE(neighbourElement)
         CLASS IS(meshCell)
           !Particle moved to new cell, repeat find procedure
-          CALL nextElement%findCell(part, self)
+          CALL neighbourElement%findCell(part, self)
 
         CLASS IS (meshEdge)
           !Particle encountered a surface, apply boundary
-          CALL nextElement%fBoundary(part%species%n)%apply(nextElement,part)
+          CALL neighbourElement%fBoundary(part%species%n)%apply(neighbourElement,part)
 
           !If particle is still inside the domain, call findCell
           IF (part%n_in) THEN 
@@ -709,7 +711,7 @@ MODULE moduleMesh
       LOGICAL:: found
       CLASS(meshCell), POINTER:: cell
       REAL(8), DIMENSION(1:3):: Xi
-      CLASS(meshElement), POINTER:: nextElement
+      CLASS(meshElement), POINTER:: neighbourElement
       INTEGER:: sp
 
       found = .FALSE.
@@ -727,11 +729,11 @@ MODULE moduleMesh
           found = .TRUE.
 
         ELSE
-          CALL cell%nextElement(Xi, nextElement)
+          CALL cell%neighbourElement(Xi, neighbourElement)
-          SELECT TYPE(nextElement)
+          SELECT TYPE(neighbourElement)
           CLASS IS(meshCell)
             !Try next element
-            cell => nextElement
+            cell => neighbourElement
 
           CLASS DEFAULT
             !Should never happend, but just in case, stops loops