diff --git a/src/makefile b/src/makefile
index 7bb277a..b398134 100644
--- a/src/makefile
+++ b/src/makefile
@@ -3,7 +3,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
 			  	$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
 			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
 			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
-			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o \
+			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o  $(OBJDIR)/moduleMath.o \
 			  	$(OBJDIR)/moduleMesh3DCart.o $(OBJDIR)/moduleMesh3DCartRead.o	\
 			  	$(OBJDIR)/moduleMesh2DCyl.o  $(OBJDIR)/moduleMesh2DCylRead.o \
 			  	$(OBJDIR)/moduleMesh2DCart.o $(OBJDIR)/moduleMesh2DCartRead.o	\
diff --git a/src/modules/makefile b/src/modules/makefile
index 28fc4ec..cea6d73 100644
--- a/src/modules/makefile
+++ b/src/modules/makefile
@@ -2,11 +2,11 @@
 OBJS = moduleCaseParam.o moduleCompTime.o moduleList.o \
 		   moduleOutput.o moduleInput.o moduleSolver.o \
 			 moduleCollisions.o moduleTable.o moduleParallel.o \
-			 moduleEM.o moduleRandom.o
+			 moduleEM.o moduleRandom.o moduleMath.o
 
 all: $(OBJS) mesh.o
 
-mesh.o: moduleCollisions.o moduleBoundary.o
+mesh.o: moduleCollisions.o moduleBoundary.o moduleMath.o
 	$(MAKE) -C mesh all
 
 moduleCollisions.o: moduleRandom.o moduleTable.o moduleSpecies.o moduleRefParam.o moduleConstParam.o moduleCollisions.f90
@@ -21,7 +21,7 @@ moduleInject.o: moduleRandom.o moduleSpecies.o moduleSolver.o moduleInject.f90
 moduleList.o: moduleSpecies.o moduleErrors.o moduleList.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
-moduleOutput.o: moduleSpecies.o moduleRefParam.o moduleOutput.f90
+moduleOutput.o: moduleMath.o moduleSpecies.o moduleRefParam.o moduleOutput.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 moduleRandom.o: moduleConstParam.o moduleRandom.f90
diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90
index 366fa84..020b971 100644
--- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90
+++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90
@@ -148,7 +148,6 @@ MODULE moduleMesh1DCart
       self%x = r1(1)
 
       self%normal = (/ 1.D0, 0.D0, 0.D0 /)
-      self%normal = self%normal/NORM2(self%normal)
 
       !Boundary index
       self%boundary => boundary(bt)
@@ -176,11 +175,11 @@ MODULE moduleMesh1DCart
 
     END FUNCTION getNodes1DCart
 
-    PURE FUNCTION intersection1DCart(self, r0, v0) RESULT(r)
+    PURE FUNCTION intersection1DCart(self, r0) RESULT(r)
       IMPLICIT NONE
 
       CLASS(meshEdge1DCart), INTENT(in):: self
-      REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
+      REAL(8), DIMENSION(1:3), INTENT(in):: r0
       REAL(8), DIMENSION(1:3):: r
 
       r = (/ self%x, 0.D0, 0.D0 /)
@@ -383,7 +382,7 @@ MODULE moduleMesh1DCart
     END FUNCTION insideSegm
 
     SUBROUTINE scatterSegm(self, part)
-      USE moduleOutput
+      USE moduleMath
       USE moduleSpecies
       IMPLICIT NONE
 
diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90
index d0f8311..08609cd 100644
--- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90
+++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90
@@ -7,6 +7,9 @@ MODULE moduleMesh1DRad
   USE moduleMeshBoundary
   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
@@ -174,11 +177,11 @@ MODULE moduleMesh1DRad
 
     END FUNCTION getNodes1DRad
 
-    PURE FUNCTION intersection1DRad(self, r0, v0) RESULT(r)
+    PURE FUNCTION intersection1DRad(self, r0) RESULT(r)
       IMPLICIT NONE
 
       CLASS(meshEdge1DRad), INTENT(in):: self
-      REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
+      REAL(8), DIMENSION(1:3), INTENT(in):: r0
       REAL(8), DIMENSION(1:3):: r
 
       r = (/ self%r, 0.D0, 0.D0 /)
@@ -317,19 +320,26 @@ MODULE moduleMesh1DRad
       REAL(8):: ke(1:2,1:2)
       REAL(8):: Xii(1:3)
       REAL(8):: dPsi(1:1, 1:2)
-      REAL(8):: invJ
+      REAL(8):: invJ(1), detJ
       REAL(8):: r, fPsi(1:2)
+      INTEGER:: l
 
       ke      = 0.D0
       Xii     = 0.D0
-      dPsi    = self%dPsi(Xii)
-      invJ    = self%invJac(Xii, dPsi)
-      fPsi    = self%fPsi(Xii)
-      r       = DOT_PRODUCT(fPsi, self%r)
-      ke(1,:) = (/ dPsi(1,1)*dPsi(1,1), dPsi(1,1)*dPsi(1,2) /) 
-      ke(2,:) = (/ dPsi(1,2)*dPsi(1,1), dPsi(1,2)*dPsi(1,2) /) 
-      ke      = 2.D0*ke*invJ
-      ke      = ke*r*PI2
+      DO l = 1, 3
+        xii(1) = corSeg(l)
+        dPsi    = self%dPsi(Xii)
+        detJ    = self%detJac(Xii, dPsi)
+        invJ    = self%invJac(Xii, dPsi)
+        fPsi    = self%fPsi(Xii)
+        r       = DOT_PRODUCT(fPsi, self%r)
+        ke      = ke + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), &
+                              RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* &
+                       r*wSeg(l)/detJ
+
+      END DO
+
+      ke      = ke*PI2
 
     END FUNCTION elemKRad
 
@@ -341,17 +351,23 @@ MODULE moduleMesh1DRad
       REAL(8), INTENT(in):: source(1:)
       REAL(8), ALLOCATABLE:: localF(:)
       REAL(8):: fPsi(1:2)
-      REAL(8):: r
-      REAL(8):: detJ
+      REAL(8):: detJ, f, r
       REAL(8):: Xii(1:3)
+      INTEGER:: l
 
-      Xii    = 0.D0
-      fPsi   = self%fPsi(Xii)
-      detJ   = self%detJac(Xii)
-      r      = DOT_PRODUCT(fPsi,self%r)
       ALLOCATE(localF(1:2))
-      localF = 2.D0*DOT_PRODUCT(fPsi, source)*detJ
-      localF = localF*r*PI2
+      localF = 0.D0
+      Xii    = 0.D0
+
+      DO l = 1, 3
+        xii(1) = corSeg(l)
+        detJ   = self%detJac(Xii)
+        fPsi   = self%fPsi(Xii)
+        r      = DOT_PRODUCT(fPsi, self%r)
+        f      = DOT_PRODUCT(fPsi, source)
+        localF = localF + f*fPsi*r*wSeg(l)*detJ
+
+      END DO
 
     END FUNCTION elemFRad
 
@@ -377,7 +393,7 @@ MODULE moduleMesh1DRad
     END FUNCTION insideRad
 
     SUBROUTINE scatterRad(self, part)
-      USE moduleOutput
+      USE moduleMath
       USE moduleSpecies
       IMPLICIT NONE
 
diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
index 40c218c..10b12fe 100644
--- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90
+++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
@@ -190,9 +190,9 @@ MODULE moduleMesh2DCart
       self%x  = (/r1(1), r2(1)/)
       self%y  = (/r1(2), r2(2)/)
       !Normal vector
-      self%normal = (/ self%y(2)-self%y(1), &
-                       self%x(2)-self%x(1), &
-                       0.D0 /)
+      self%normal = (/ -(self%y(2)-self%y(1)), &
+                         self%x(2)-self%x(1) , &
+                         0.D0 /)
       self%normal = self%normal/NORM2(self%normal)
 
       !Boundary index
@@ -243,20 +243,21 @@ MODULE moduleMesh2DCart
 
     END FUNCTION getNodes2DCart
 
-    PURE FUNCTION intersection2DCartEdge(self, r0, v0) RESULT(r)
+    PURE FUNCTION intersection2DCartEdge(self, r0) RESULT(r)
       IMPLICIT NONE
 
       CLASS(meshEdge2DCart), INTENT(in):: self
-      REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
+      REAL(8), DIMENSION(1:3), INTENT(in):: r0
       REAL(8), DIMENSION(1:3):: r
-      REAL(8), DIMENSION(1:3):: rS !base point of surface
-      REAL(8):: d
+      REAL(8), DIMENSION(1:3):: edge0, edgeV
+      REAL(8):: tI
 
-      rS = (/ self%x(1), self%y(1), 0.D0 /)
+      edge0 = (/self%x(1), self%y(1), 0.D0 /)
+      edgeV = (/self%x(2), self%y(2), 0.D0 /) - edge0
 
-      d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal)
+      tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV)
 
-      r = r0 + v0*d
+      r = edge0 + tI*edgeV
 
     END FUNCTION intersection2DCartEdge
 
@@ -496,7 +497,7 @@ MODULE moduleMesh2DCart
 
     !Scatter properties of particle into element nodes
     SUBROUTINE scatterQuad(self, part)
-      USE moduleOutput
+      USE moduleMath
       USE moduleSpecies
       IMPLICIT NONE
 
@@ -855,7 +856,7 @@ MODULE moduleMesh2DCart
 
     !Scatter properties of particles into element
     SUBROUTINE scatterTria(self, part)
-      USE moduleOutput
+      USE moduleMath
       USE moduleSpecies
       IMPLICIT NONE
 
diff --git a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90 b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90
index b0f2ee3..96664df 100644
--- a/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90
+++ b/src/modules/mesh/2DCart/moduleMesh2DCartRead.f90
@@ -437,15 +437,15 @@ MODULE moduleMesh2DCartRead
         IF     (elemA%n1%n == elemB%n1%n .AND. &
                 elemA%n2%n == elemB%n2%n) THEN
           elemA%e1 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
                 elemA%n2%n == elemB%n1%n) THEN
           elemA%e1 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -456,15 +456,15 @@ MODULE moduleMesh2DCartRead
         IF     (elemA%n2%n == elemB%n1%n .AND. &
                 elemA%n3%n == elemB%n2%n) THEN
           elemA%e2 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
                 elemA%n3%n == elemB%n1%n) THEN
           elemA%e2 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -475,15 +475,15 @@ MODULE moduleMesh2DCartRead
         IF     (elemA%n3%n == elemB%n1%n .AND. &
                 elemA%n4%n == elemB%n2%n) THEN
           elemA%e3 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
                 elemA%n4%n == elemB%n1%n) THEN
           elemA%e3 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -494,15 +494,15 @@ MODULE moduleMesh2DCartRead
         IF     (elemA%n4%n == elemB%n1%n .AND. &
                 elemA%n1%n == elemB%n2%n) THEN
           elemA%e4 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
                 elemA%n1%n == elemB%n1%n) THEN
           elemA%e4 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -521,16 +521,16 @@ MODULE moduleMesh2DCartRead
         IF     (elemA%n1%n == elemB%n1%n .AND. &
                 elemA%n2%n == elemB%n2%n) THEN
           elemA%e1 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
                 elemA%n2%n == elemB%n1%n) THEN
           elemA%e1 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
 
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
+          
         END IF
 
       END IF
@@ -540,15 +540,15 @@ MODULE moduleMesh2DCartRead
         IF     (elemA%n2%n == elemB%n1%n .AND. &
                 elemA%n3%n == elemB%n2%n) THEN
           elemA%e2 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
                 elemA%n3%n == elemB%n1%n) THEN
           elemA%e2 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -559,15 +559,15 @@ MODULE moduleMesh2DCartRead
         IF     (elemA%n3%n == elemB%n1%n .AND. &
                 elemA%n1%n == elemB%n2%n) THEN
           elemA%e3 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
                 elemA%n1%n == elemB%n1%n) THEN
           elemA%e3 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
index e0f3a27..d48af83 100644
--- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
+++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
@@ -191,9 +191,9 @@ MODULE moduleMesh2DCyl
       self%z  = (/r1(1), r2(1)/)
       self%r  = (/r1(2), r2(2)/)
       !Normal vector
-      self%normal = (/ self%r(2)-self%r(1), &
-                       self%z(2)-self%z(1), &
-                       0.D0 /)
+      self%normal = (/ -(self%r(2)-self%r(1)), &
+                         self%z(2)-self%z(1) , &
+                         0.D0 /)
       self%normal = self%normal/NORM2(self%normal)
       !Boundary index
       self%boundary => boundary(bt)
@@ -221,20 +221,22 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION getNodes2DCyl
 
-    PURE FUNCTION intersection2DCylEdge(self, r0, v0) RESULT(r)
+    PURE FUNCTION intersection2DCylEdge(self, r0) RESULT(r)
+      USE moduleMath
       IMPLICIT NONE
 
       CLASS(meshEdge2DCyl), INTENT(in):: self
-      REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
+      REAL(8), DIMENSION(1:3), INTENT(in):: r0
       REAL(8), DIMENSION(1:3):: r
-      REAL(8), DIMENSION(1:3):: rS !base point of surface
-      REAL(8):: d
+      REAL(8), DIMENSION(1:3):: edge0, edgeV
+      REAL(8):: tI
 
-      rS = (/ self%z(1), self%r(1), 0.D0 /)
+      edge0 = (/self%z(1), self%r(1), 0.D0 /)
+      edgeV = (/self%z(2), self%r(2), 0.D0 /) - edge0
 
-      d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal)
+      tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV)
 
-      r = r0 + v0*d
+      r = edge0 + tI*edgeV
 
     END FUNCTION intersection2DCylEdge
 
@@ -507,7 +509,7 @@ MODULE moduleMesh2DCyl
 
     !Scatter properties of particle into element nodes
     SUBROUTINE scatterQuad(self, part)
-      USE moduleOutput
+      USE moduleMath
       USE moduleSpecies
       IMPLICIT NONE
 
@@ -874,7 +876,7 @@ MODULE moduleMesh2DCyl
 
     !Scatter properties of particles into element
     SUBROUTINE scatterTria(self, part)
-      USE moduleOutput
+      USE moduleMath
       USE moduleSpecies
       IMPLICIT NONE
 
diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90
index 379a0fb..cbe3faf 100644
--- a/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90
+++ b/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90
@@ -217,6 +217,29 @@ MODULE moduleMesh2DCylRead
 
     END SUBROUTINE connectedVolEdge
 
+    PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident)
+      IMPLICIT NONE
+
+      INTEGER, DIMENSION(1:2), INTENT(in):: nodesA, nodesB
+      LOGICAL:: coincident
+      INTEGER:: i
+
+      coincident = .FALSE.
+      DO i = 1, 2
+        IF (ANY(nodesA(i) == nodesB)) THEN
+          coincident = .TRUE.
+
+        ELSE
+          coincident = .FALSE.
+          EXIT
+
+        END IF
+
+      END DO
+
+    END FUNCTION coincidentNodes
+
+
     SUBROUTINE connectedQuadQuad(elemA, elemB)
       IMPLICIT NONE
 
@@ -437,15 +460,15 @@ MODULE moduleMesh2DCylRead
         IF     (elemA%n1%n == elemB%n1%n .AND. &
                 elemA%n2%n == elemB%n2%n) THEN
           elemA%e1 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
                 elemA%n2%n == elemB%n1%n) THEN
           elemA%e1 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -456,15 +479,15 @@ MODULE moduleMesh2DCylRead
         IF     (elemA%n2%n == elemB%n1%n .AND. &
                 elemA%n3%n == elemB%n2%n) THEN
           elemA%e2 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
                 elemA%n3%n == elemB%n1%n) THEN
           elemA%e2 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -475,15 +498,15 @@ MODULE moduleMesh2DCylRead
         IF     (elemA%n3%n == elemB%n1%n .AND. &
                 elemA%n4%n == elemB%n2%n) THEN
           elemA%e3 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
                 elemA%n4%n == elemB%n1%n) THEN
           elemA%e3 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -494,12 +517,15 @@ MODULE moduleMesh2DCylRead
         IF     (elemA%n4%n == elemB%n1%n .AND. &
                 elemA%n1%n == elemB%n2%n) THEN
           elemA%e4 => elemB
-          elemB%e2 => elemA
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
                 elemA%n1%n == elemB%n1%n) THEN
           elemA%e4 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -518,15 +544,15 @@ MODULE moduleMesh2DCylRead
         IF     (elemA%n1%n == elemB%n1%n .AND. &
                 elemA%n2%n == elemB%n2%n) THEN
           elemA%e1 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
                 elemA%n2%n == elemB%n1%n) THEN
           elemA%e1 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -537,15 +563,15 @@ MODULE moduleMesh2DCylRead
         IF     (elemA%n2%n == elemB%n1%n .AND. &
                 elemA%n3%n == elemB%n2%n) THEN
           elemA%e2 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
                 elemA%n3%n == elemB%n1%n) THEN
           elemA%e2 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
@@ -556,15 +582,15 @@ MODULE moduleMesh2DCylRead
         IF     (elemA%n3%n == elemB%n1%n .AND. &
                 elemA%n1%n == elemB%n2%n) THEN
           elemA%e3 => elemB
-          elemB%e2 => elemA
-
-          !Revers the normal to point inside the domain
-          elemB%normal = - elemB%normal
+          elemB%e1 => elemA
 
         ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
                 elemA%n1%n == elemB%n1%n) THEN
           elemA%e3 => elemB
-          elemB%e1 => elemA
+          elemB%e2 => elemA
+
+          !Revers the normal to point inside the domain
+          elemB%normal = - elemB%normal
 
         END IF
 
diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90
index 315cb42..ec7f2b1 100644
--- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90
+++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90
@@ -134,6 +134,7 @@ MODULE moduleMesh3DCart
       USE moduleSpecies
       USE moduleBoundary
       USE moduleErrors
+      USE moduleMath
       IMPLICIT NONE
 
       CLASS(meshEdge3DCartTria), INTENT(out):: self
@@ -142,6 +143,7 @@ MODULE moduleMesh3DCart
       INTEGER, INTENT(in):: bt
       INTEGER, INTENT(in):: physicalSurface
       REAL(8), DIMENSION(1:3):: r1, r2, r3
+      REAL(8), DIMENSION(1:3):: vec1, vec2
       INTEGER:: s
 
       self%n = n
@@ -156,10 +158,14 @@ MODULE moduleMesh3DCart
       self%y  = (/r1(2), r2(2), r3(2)/)
       self%z  = (/r1(3), r2(3), r3(3)/)
       !Normal vector
-      self%normal = (/  (self%y(2)-self%y(1))*(self%z(3)-self%z(1)) - (self%z(2)-self%z(1))*(self%y(3)-self%y(1)), &
-                        (self%x(2)-self%x(1))*(self%z(3)-self%z(1)) - (self%z(2)-self%z(1))*(self%x(3)-self%x(1)), &
-                        (self%x(2)-self%x(1))*(self%y(3)-self%y(1)) - (self%z(2)-self%z(1))*(self%y(3)-self%y(1)) /)
-      self%normal = self%normal/NORM2(self%normal)
+      vec1 = (/ self%x(2) - self%x(1), &
+                self%y(2) - self%y(1), &
+                self%z(2) - self%z(1) /)
+      vec2 = (/ self%x(3) - self%x(1), &
+                self%y(3) - self%y(1), &
+                self%z(3) - self%z(1) /)
+      self%normal = crossProduct(vec1, vec2)
+      self%normal = normalize(self%normal)
 
       !Boundary index
       self%boundary => boundary(bt)
@@ -187,20 +193,21 @@ MODULE moduleMesh3DCart
 
     END FUNCTION getNodes3DCartTria
 
-    PURE FUNCTION intersection3DCartTria(self, r0, v0) RESULT(r)
+    PURE FUNCTION intersection3DCartTria(self, r0) RESULT(r)
       IMPLICIT NONE
 
       CLASS(meshEdge3DCartTria), INTENT(in):: self
-      REAL(8), DIMENSION(1:3), INTENT(in):: r0, v0
+      REAL(8), INTENT(in):: r0(1:3)
       REAL(8), DIMENSION(1:3):: r
-      REAL(8), DIMENSION(1:3):: rS !base point of surface
-      REAL(8):: d
+      REAL(8), DIMENSION(1:3):: edge0, edgeV
+      REAL(8):: tI
 
-      rS = (/ self%x(1), self%y(1), self%z(1) /)
+      edge0 = (/self%x(1), self%y(1), self%z(1) /)
+      edgeV = (/self%x(2), self%y(2), self%z(2) /) - edge0
 
-      d = DOT_PRODUCT((rS - r0), self%normal)/DOT_PRODUCT(v0, self%normal)
+      tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV)
 
-      r = r0 + v0*d
+      r = edge0 + tI*edgeV
 
     END FUNCTION intersection3DCartTria
 
@@ -476,7 +483,7 @@ MODULE moduleMesh3DCart
     END FUNCTION insideTetra
 
     SUBROUTINE scatterTetra(self, part)
-      USE moduleOutput
+      USE moduleMath
       USE moduleSpecies
       IMPLICIT NONE
 
@@ -577,7 +584,7 @@ MODULE moduleMesh3DCart
       INTEGER:: nextInt
 
       !TODO: Review when connectivity
-      xiArray = (/ xi(3), xi(2), 1.D0 - xi(1) - xi(2) - xi(3), xi(1) /)
+      xiArray = (/ xi(3), 1.D0 - xi(1) - xi(2) - xi(3), xi(2), xi(1) /)
       nextInt = MINLOC(xiArray, 1)
       NULLIFY(nextElement)
       SELECT CASE(nextInt)
@@ -647,9 +654,11 @@ MODULE moduleMesh3DCart
       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(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)
 
     END FUNCTION invJ3DCart
 
diff --git a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90 b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90
index c6b021a..01aef8e 100644
--- a/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90
+++ b/src/modules/mesh/3DCart/moduleMesh3DCartRead.f90
@@ -269,25 +269,25 @@ MODULE moduleMesh3DCartRead
 
       !Check surface 2
       IF (.NOT. ASSOCIATED(elemA%e2)) THEN
-        IF     (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
+        IF     (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
                                 (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
 
           elemA%e2 => elemB
           elemB%e1 => elemA
 
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
+        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
                                 (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
 
           elemA%e2 => elemB
           elemB%e2 => elemA
 
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
+        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
                                 (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
 
           elemA%e2 => elemB
           elemB%e3 => elemA
 
-        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
+        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
                                 (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
 
           elemA%e2 => elemB
@@ -299,25 +299,25 @@ MODULE moduleMesh3DCartRead
 
       !Check surface 3
       IF (.NOT. ASSOCIATED(elemA%e3)) THEN
-        IF     (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+        IF     (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
                                 (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
 
           elemA%e3 => elemB
           elemB%e1 => elemA
 
-        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
                                 (/elemB%n2%n, elemB%n3%n, elemB%n4%n/))) THEN
 
           elemA%e3 => elemB
           elemB%e2 => elemA
 
-        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
                                 (/elemB%n1%n, elemB%n2%n, elemB%n4%n/))) THEN
 
           elemA%e3 => elemB
           elemB%e3 => elemA
 
-        ELSEIF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+        ELSEIF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elemA%n4%n/), &
                                 (/elemB%n1%n, elemB%n3%n, elemB%n4%n/))) THEN
 
           elemA%e3 => elemB
@@ -328,7 +328,7 @@ MODULE moduleMesh3DCartRead
       END IF
 
       !Check surface 4
-      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
         IF     (coincidentNodes((/elemA%n1%n, elemA%n3%n, elemA%n4%n/), &
                                 (/elemB%n1%n, elemB%n2%n, elemB%n3%n/))) THEN
 
@@ -360,77 +360,78 @@ MODULE moduleMesh3DCartRead
     END SUBROUTINE connectedTetraTetra
 
     SUBROUTINE connectedTetraEdge(elemA, elemB)
+      USE moduleMath
       IMPLICIT NONE
 
       CLASS(meshVol3DCartTetra), INTENT(inout), TARGET:: elemA
       CLASS(meshEdge3DCartTria), INTENT(inout), TARGET:: elemB
+      INTEGER:: nodesEdge(1:3)
+      REAL(8), DIMENSION(1:3):: vec1, vec2
+      REAL(8):: normVol(1:3)
+
+      nodesEdge = (/ elemB%n1%n, elemB%n2%n, elemB%n3%n /)
 
       !Check surface 1
       IF (.NOT. ASSOCIATED(elemA%e1)) THEN
-        IF     ((elemA%n1%n == elemB%n1%n .AND. &
-                 elemA%n2%n == elemB%n2%n .AND. &
-                 elemA%n3%n == elemB%n3%n) .OR. &
-                (elemA%n1%n == elemB%n3%n .AND. &
-                 elemA%n2%n == elemB%n1%n .AND. &
-                 elemA%n3%n == elemB%n2%n) .OR. &
-                (elemA%n1%n == elemB%n2%n .AND. &
-                 elemA%n2%n == elemB%n3%n .AND. &
-                 elemA%n3%n == elemB%n1%n)) THEN
+        IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n3%n/), &
+                              nodesEdge)) THEN
 
-          elemA%e1 => elemB
-          elemB%e1 => elemA
+          vec1 = (/ elemA%x(2) - elemA%x(1), &
+                    elemA%y(2) - elemA%y(1), &
+                    elemA%z(2) - elemA%z(1) /)
+          vec2 = (/ elemA%x(3) - elemA%x(1), &
+                    elemA%y(3) - elemA%y(1), &
+                    elemA%z(3) - elemA%z(1) /)
+          normVol = crossProduct(vec1, vec2)
+          normVol = normalize(normVol)
 
-        ELSEIF ((elemA%n1%n == elemB%n3%n .AND. &
-                 elemA%n2%n == elemB%n2%n .AND. &
-                 elemA%n3%n == elemB%n1%n) .OR. &
-                (elemA%n1%n == elemB%n1%n .AND. &
-                 elemA%n2%n == elemB%n3%n .AND. &
-                 elemA%n3%n == elemB%n2%n) .OR. &
-                (elemA%n1%n == elemB%n2%n .AND. &
-                 elemA%n2%n == elemB%n1%n .AND. &
-                 elemA%n3%n == elemB%n3%n)) THEN
+          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
 
-          elemA%e1 => elemB
-          elemB%e2 => elemA
+            elemA%e1 => elemB
+            elemB%e1 => elemA
 
-          !Revers the normal to point inside the domain
-          elemB%normal = -elemB%normal
+          ELSE
+
+            elemA%e1 => elemB
+            elemB%e2 => elemA
+
+            !Revers the normal to point inside the domain
+            elemB%normal = -elemB%normal
+
+          END IF
 
         END IF
 
-
       END IF
 
       !Check surface 2
       IF (.NOT. ASSOCIATED(elemA%e2)) THEN
-        IF     ((elemA%n1%n == elemB%n1%n .AND. &
-                 elemA%n2%n == elemB%n2%n .AND. &
-                 elemA%n4%n == elemB%n3%n) .OR. &
-                (elemA%n1%n == elemB%n3%n .AND. &
-                 elemA%n2%n == elemB%n1%n .AND. &
-                 elemA%n4%n == elemB%n2%n) .OR. &
-                (elemA%n1%n == elemB%n2%n .AND. &
-                 elemA%n2%n == elemB%n3%n .AND. &
-                 elemA%n4%n == elemB%n1%n)) THEN
+        IF (coincidentNodes((/elemA%n2%n, elemA%n3%n, elemA%n4%n/), &
+                              nodesEdge)) THEN
 
-          elemA%e2 => elemB
-          elemB%e1 => elemA
+          vec1 = (/ elemA%x(3) - elemA%x(2), &
+                    elemA%y(3) - elemA%y(2), &
+                    elemA%z(3) - elemA%z(2) /)
+          vec2 = (/ elemA%x(4) - elemA%x(2), &
+                    elemA%y(4) - elemA%y(2), &
+                    elemA%z(4) - elemA%z(2) /)
+          normVol = crossProduct(vec1, vec2)
+          normVol = normalize(normVol)
 
-        ELSEIF ((elemA%n1%n == elemB%n3%n .AND. &
-                 elemA%n2%n == elemB%n2%n .AND. &
-                 elemA%n4%n == elemB%n1%n) .OR. &
-                (elemA%n1%n == elemB%n1%n .AND. &
-                 elemA%n2%n == elemB%n3%n .AND. &
-                 elemA%n4%n == elemB%n2%n) .OR. &
-                (elemA%n1%n == elemB%n2%n .AND. &
-                 elemA%n2%n == elemB%n1%n .AND. &
-                 elemA%n4%n == elemB%n3%n)) THEN
+          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
 
-          elemA%e2 => elemB
-          elemB%e2 => elemA
+            elemA%e2 => elemB
+            elemB%e1 => elemA
 
-          !Revers the normal to point inside the domain
-          elemB%normal = -elemB%normal
+          ELSE
+
+            elemA%e2 => elemB
+            elemB%e2 => elemA
+
+            !Revers the normal to point inside the domain
+            elemB%normal = -elemB%normal
+
+          END IF
 
         END IF
 
@@ -438,66 +439,66 @@ MODULE moduleMesh3DCartRead
 
       !Check surface 3
       IF (.NOT. ASSOCIATED(elemA%e3)) THEN
-        IF     ((elemA%n2%n == elemB%n1%n .AND. &
-                 elemA%n3%n == elemB%n2%n .AND. &
-                 elemA%n4%n == elemB%n3%n) .OR. &
-                (elemA%n2%n == elemB%n3%n .AND. &
-                 elemA%n3%n == elemB%n1%n .AND. &
-                 elemA%n4%n == elemB%n2%n) .OR. &
-                (elemA%n2%n == elemB%n2%n .AND. &
-                 elemA%n3%n == elemB%n3%n .AND. &
-                 elemA%n4%n == elemB%n1%n)) THEN
+        IF (coincidentNodes((/elemA%n1%n, elemA%n2%n, elema%n4%n/), &
+                              nodesEdge)) THEN
 
-          elemA%e3 => elemB
-          elemB%e1 => elemA
+          vec1 = (/ elemA%x(2) - elemA%x(1), &
+                    elemA%y(2) - elemA%y(1), &
+                    elemA%z(2) - elemA%z(1) /)
+          vec2 = (/ elemA%x(4) - elemA%x(1), &
+                    elemA%y(4) - elemA%y(1), &
+                    elemA%z(4) - elemA%z(1) /)
+          normVol = crossProduct(vec1, vec2)
+          normVol = normalize(normVol)
 
-        ELSEIF ((elemA%n2%n == elemB%n3%n .AND. &
-                 elemA%n3%n == elemB%n2%n .AND. &
-                 elemA%n4%n == elemB%n1%n) .OR. &
-                (elemA%n2%n == elemB%n1%n .AND. &
-                 elemA%n3%n == elemB%n3%n .AND. &
-                 elemA%n4%n == elemB%n2%n) .OR. &
-                (elemA%n2%n == elemB%n2%n .AND. &
-                 elemA%n3%n == elemB%n1%n .AND. &
-                 elemA%n4%n == elemB%n3%n)) THEN
+          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
 
-          elemA%e3 => elemB
-          elemB%e2 => elemA
+            elemA%e3 => elemB
+            elemB%e1 => elemA
 
-          !Revers the normal to point inside the domain
-          elemB%normal = -elemB%normal
+          ELSE
+
+            elemA%e3 => elemB
+            elemB%e2 => elemA
+
+            !Revers the normal to point inside the domain
+            elemB%normal = -elemB%normal
+
+          END IF
 
         END IF
 
       END IF
 
       !Check surface 4
-      IF (.NOT. ASSOCIATED(elemA%e3)) THEN
-        IF     ((elemA%n1%n == elemB%n1%n .AND. &
-                 elemA%n3%n == elemB%n2%n .AND. &
-                 elemA%n4%n == elemB%n3%n) .OR. &
-                (elemA%n1%n == elemB%n3%n .AND. &
-                 elemA%n3%n == elemB%n1%n .AND. &
-                 elemA%n4%n == elemB%n2%n) .OR. &
-                (elemA%n1%n == elemB%n2%n .AND. &
-                 elemA%n3%n == elemB%n3%n .AND. &
-                 elemA%n4%n == elemB%n1%n)) THEN
+      IF (.NOT. ASSOCIATED(elemA%e4)) THEN
+        IF (coincidentNodes((/elemA%n1%n, elemA%n3%n, elema%n4%n/), &
+                              nodesEdge)) THEN
 
-          elemA%e4 => elemB
-          elemB%e1 => elemA
+          vec1 = (/ elemA%x(3) - elemA%x(1), &
+                    elemA%y(3) - elemA%y(1), &
+                    elemA%z(3) - elemA%z(1) /)
+          vec2 = (/ elemA%x(4) - elemA%x(1), &
+                    elemA%y(4) - elemA%y(1), &
+                    elemA%z(4) - elemA%z(1) /)
+          normVol = crossProduct(vec1, vec2)
+          normVol = normalize(normVol)
 
-        ELSEIF ((elemA%n1%n == elemB%n3%n .AND. &
-                 elemA%n3%n == elemB%n2%n .AND. &
-                 elemA%n4%n == elemB%n1%n) .OR. &
-                (elemA%n1%n == elemB%n1%n .AND. &
-                 elemA%n3%n == elemB%n3%n .AND. &
-                 elemA%n4%n == elemB%n2%n) .OR. &
-                (elemA%n1%n == elemB%n2%n .AND. &
-                 elemA%n3%n == elemB%n1%n .AND. &
-                 elemA%n4%n == elemB%n3%n)) THEN
+          IF (DOT_PRODUCT(elemB%normal, normVol) == -1.D0) THEN
 
-          elemA%e4 => elemB
-          elemB%e2 => elemA
+            elemA%e4 => elemB
+            elemB%e1 => elemA
+
+
+          ELSE
+
+            elemA%e4 => elemB
+            elemB%e2 => elemA
+
+            !Revers the normal to point inside the domain
+            elemB%normal = -elemB%normal
+
+          END IF
 
         END IF
 
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index f5f7b47..137b256 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -97,10 +97,10 @@ MODULE moduleMesh
     END FUNCTION getNodesEdge_interface
 
     !Returns the intersecction between an edge and a line defined by point r0 and vector v0
-    PURE FUNCTION intersectionEdge_interface(self, r0, v0) RESULT(r)
+    PURE FUNCTION intersectionEdge_interface(self, r0) RESULT(r)
       IMPORT:: meshEdge
       CLASS(meshEdge), INTENT(in):: self
-      REAL(8), INTENT(in), DIMENSION(1:3):: r0, v0
+      REAL(8), INTENT(in), DIMENSION(1:3):: r0
       REAL(8):: r(1:3)
 
     END FUNCTION intersectionEdge_interface
diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90
index 2252b84..7f5fb88 100644
--- a/src/modules/mesh/moduleMeshBoundary.f90
+++ b/src/modules/mesh/moduleMeshBoundary.f90
@@ -13,30 +13,26 @@ MODULE moduleMeshBoundary
       !rp  = intersection between particle and edge
       !rpp = final position of particle
       !vpp = final velocity of particle
-      REAL(8), DIMENSION(1:3):: rp, rpp, vpp
+      REAL(8), DIMENSION(1:3):: rp, vpp
+      REAL(8):: tI
       REAL(8):: taup !time step for reflecting process
 
       !Reflect particle velocity
       vpp = part%v - 2.D0*DOT_PRODUCT(part%v, edge%normal)*edge%normal
+      part%v = vpp
 
-      !Computes the intersection between particle and surface
-      rp = edge%intersection(part%r, part%v)
+      rp = edge%intersection(part%r)
 
-      !Computes the reflection time step
-      taup = NORM2(part%r - rp)*tau(part%sp)
+      part%r = 2.D0*(rp - part%r) + part%r
 
-      !New position of particle
-      rpp = rp + vpp*taup
-
-      !assign new parameters to particle
-      part%r    = rpp
-      part%v    = vpp
+      !particle is assumed to be inside
       part%n_in = .TRUE.
 
     END SUBROUTINE reflection
 
     !Absoption in a surface
     SUBROUTINE absorption(edge, part)
+      USE moduleCaseParam
       USE moduleSpecies
       IMPLICIT NONE
 
@@ -45,7 +41,7 @@ MODULE moduleMeshBoundary
       REAL(8):: rpp(1:3) !Position of particle projected to the edge
       REAL(8):: d !Distance from particle to edge
 
-      rpp = edge%intersection(part%r, part%v)
+      rpp = edge%intersection(part%r)
 
       d = NORM2(rpp - part%r)
 
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index df9df48..9d50abc 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -200,6 +200,7 @@ MODULE moduleInject
       USE moduleSolver
       USE moduleMesh
       USE moduleRandom
+      USE moduleErrors
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self
@@ -241,7 +242,7 @@ MODULE moduleInject
           partInj(n)%vol = randomEdge%e2%n
 
         ELSE
-          PRINT *, "ERROR NO VOL ASSOCIATED TO EDGE"
+          CALL criticalError("No Volume associated to edge", 'addParticles') 
 
         END IF
 
diff --git a/src/modules/moduleMath.f90 b/src/modules/moduleMath.f90
new file mode 100644
index 0000000..213afa5
--- /dev/null
+++ b/src/modules/moduleMath.f90
@@ -0,0 +1,42 @@
+MODULE moduleMath
+  IMPLICIT NONE
+
+  CONTAINS
+    !Outer product of two tensors
+    PURE FUNCTION outerProduct(a,b) RESULT(s)
+      IMPLICIT NONE
+
+      REAL(8), DIMENSION(1:3), INTENT(in):: a, b
+      REAL(8):: s(1:3,1:3)
+
+      s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3)
+
+    END FUNCTION outerProduct
+
+    !Cross product of two 3D vectors
+    PURE FUNCTION crossProduct(a, b) RESULT(c)
+      IMPLICIT NONE
+
+      REAL(8), DIMENSION(1:3), INTENT(in):: a, b
+      REAL(8), DIMENSION(1:3):: c
+
+      c = 0.D0
+
+      c(1) =   a(2)*b(3) - a(3)*b(2)
+      c(2) = -(a(1)*b(3) - a(3)*b(1))
+      c(3) =   a(1)*b(2) - a(2)*b(1)
+
+    END FUNCTION crossProduct
+
+    !Normalizes a 3D vector
+    PURE FUNCTION normalize(a) RESULT(b)
+      IMPLICIT NONE
+
+      REAL(8), DIMENSION(1:3), INTENT(in):: a
+      REAL(8), DIMENSION(1:3):: b
+      
+      b = a / NORM2(a)
+
+    END FUNCTION normalize
+
+END MODULE moduleMath
diff --git a/src/modules/moduleOutput.f90 b/src/modules/moduleOutput.f90
index b4b8317..963302d 100644
--- a/src/modules/moduleOutput.f90
+++ b/src/modules/moduleOutput.f90
@@ -31,16 +31,6 @@ MODULE moduleOutput
   LOGICAL:: emOutput   = .FALSE.
 
   CONTAINS
-    FUNCTION outerProduct(a,b) RESULT(s)
-      IMPLICIT NONE
-
-      REAL(8), DIMENSION(1:3):: a, b
-      REAL(8):: s(1:3,1:3)
-
-      s = SPREAD(a, dim = 2, ncopies = 3)*SPREAD(b, dim = 1, ncopies = 3)
-
-    END FUNCTION outerProduct
-
     FUNCTION tensorTrace(a) RESULT(t)
       IMPLICIT NONE
 
@@ -56,6 +46,7 @@ MODULE moduleOutput
       USE moduleConstParam
       USE moduleRefParam
       USE moduleSpecies
+      USE moduleMath
       IMPLICIT NONE
 
       TYPE(outputNode),   INTENT(in)::  rawValues