diff --git a/src/modules/moduleMeshCyl.f95 b/src/modules/moduleMeshCyl.f95
index 7f3a4f3..d5c464f 100644
--- a/src/modules/moduleMeshCyl.f95
+++ b/src/modules/moduleMeshCyl.f95
@@ -108,6 +108,8 @@ MODULE moduleMeshCyl
     !Connectivity to adjacent elements
     CLASS(*), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()
     REAL(8):: arNodes(1:3) = 0.D0
+    !Derivatives in z,r real space
+    REAL(8), DIMENSION(1:3):: dPsiZ, dPsiR
 
     CONTAINS
       PROCEDURE, PASS::   init        => initVolTriaCyl
@@ -618,39 +620,16 @@ MODULE moduleMeshCyl
       CLASS(meshVolCylTria), INTENT(out):: self
       INTEGER, INTENT(in):: n
       INTEGER, INTENT(in):: p(:)
-      REAL(8):: ldCorner(1:3) !left down (ld) corner of square enclousing the triangle
-      INTEGER:: minP
-      INTEGER:: pOrder(1:3)
       REAL(8), DIMENSION(1:3):: r1, r2, r3
+      REAL(8):: A
 
       !Assign node index
       self%n = n
 
-      !Find the first node
-      r1 = mesh%nodes(p(1))%obj%getCoordinates()
-      r2 = mesh%nodes(p(2))%obj%getCoordinates()
-      r3 = mesh%nodes(p(3))%obj%getCoordinates()
-      ldCorner =  (/ MINVAL((/r1(1), r2(1), r3(1)/), 1), &
-                     MINVAL((/r1(2), r2(2), r3(2)/), 1), &
-                     0.D0 /)
-      minP = MINLOC((/ NORM2(r1-ldCorner), NORM2(r2-ldCorner), NORM2(r3-ldCorner) /), 1)
-      pOrder(1) = p(minP)
-      SELECT CASE(minP)
-      CASE(1)
-        pOrder(2) = p(2)
-        pOrder(3) = p(3)
-      CASE(2)
-        pOrder(2) = p(3)
-        pOrder(3) = p(1)
-      CASE(3)
-        pOrder(2) = p(1)
-        pOrder(3) = p(2)
-      END SELECT
-
       !Assign nodes to element
-      self%n1 => mesh%nodes(pOrder(1))%obj
-      self%n2 => mesh%nodes(pOrder(2))%obj
-      self%n3 => mesh%nodes(pOrder(3))%obj
+      self%n1 => mesh%nodes(p(1))%obj
+      self%n2 => mesh%nodes(p(2))%obj
+      self%n3 => mesh%nodes(p(3))%obj
       !Get element coordinates
       r1 = self%n1%getCoordinates()
       r2 = self%n2%getCoordinates()
@@ -663,6 +642,20 @@ MODULE moduleMeshCyl
       self%n2%v = self%n2%v + self%arNodes(2)
       self%n3%v = self%n3%v + self%arNodes(3)
 
+      !Derivatives in z/r for shape functions (node independent)
+      !TODO: This is used because invJ.dPsi does not produce the right Electric field
+      A = self%z(2)*self%r(3) - self%z(3)*self%r(2) + &
+          self%z(3)*self%r(1) - self%z(1)*self%r(3) + &
+          self%z(1)*self%r(2) - self%z(2)*self%r(1)
+
+      self%dPsiZ = (/ self%r(2)-self%r(3), &
+                      self%r(3)-self%r(1), &
+                      self%r(1)-self%r(2) /)/A
+
+      self%dPsiR = (/ self%z(3)-self%z(2), &
+                      self%z(1)-self%z(3), &
+                      self%z(2)-self%z(1) /)/A
+
       self%sigmaVrelMax = sigma_ref/L_ref**2
 
       CALL OMP_INIT_LOCK(self%lock)
@@ -882,9 +875,6 @@ MODULE moduleMeshCyl
 
       CLASS(meshVolCylTria), INTENT(in):: self
       REAL(8), INTENT(in):: xi(1:3)
-      REAL(8):: dPsi(1:2,1:3)
-      REAL(8):: dPsiR(1:2,1:3)
-      REAL(8):: invJ(1:2,1:2), detJ
       REAL(8):: phi(1:3)
       REAL(8):: EF(1:3)
 
@@ -892,12 +882,8 @@ MODULE moduleMeshCyl
               self%n2%emData%phi, &
               self%n3%emData%phi/)
 
-      dPsi  =  self%dPsi(xi)
-      detJ  =  self%detJac(xi,dPsi)
-      invJ  =  self%invJac(xi,dPsi)
-      dPsiR =  MATMUL(invJ, dPsi)/detJ
-      EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi)
-      EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi)
+      EF(1) = -DOT_PRODUCT(self%dPsiZ(:), phi)
+      EF(2) = -DOT_PRODUCT(self%dPsiR(:), phi)
       EF(3) =  0.D0
 
     END FUNCTION gatherEFTria
@@ -926,7 +912,7 @@ MODULE moduleMeshCyl
       REAL(8):: dPsi(1:2,1:3)
 
       !Direct method to convert coordinates
-      xi      = (/ 0.D0, 0.D0, 0.D0 /) !Irrelevant, required for input
+      xi      = 0.D0 !Irrelevant, required for input
       deltaR  = (/ r(1) - self%z(1), r(2) - self%r(1) /)
       dPsi    = self%dPsi(xi)
       invJ    = self%invJac(xi, dPsi)