diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
index 426846d..3f88902 100644
--- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
+++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
@@ -144,6 +144,7 @@ MODULE moduleMesh2DCyl
       USE moduleSpecies
       USE moduleBoundary
       USE moduleErrors
+      USE moduleConstParam, ONLY: PI
       IMPLICIT NONE
 
       CLASS(meshEdge2DCyl), INTENT(out):: self
@@ -163,7 +164,15 @@ MODULE moduleMesh2DCyl
       r2 = self%n2%getCoordinates()
       self%z  = (/r1(1), r2(1)/)
       self%r  = (/r1(2), r2(2)/)
-      self%weight = DABS(self%r(2)**2 - self%r(1)**2)
+      !Edge surface
+      IF (self%z(2) /= self%z(1)) THEN
+        self%surface = ABS(self%r(2) + self%r(1))*ABS(self%z(2) - self%z(1))
+
+      ELSE
+        self%surface = ABS(self%r(2)**2 - self%r(1)**2)
+
+      END IF
+      self%surface = self%surface * PI
       !Normal vector
       self%normal = (/ -(self%r(2)-self%r(1)), &
                          self%z(2)-self%z(1) , &
@@ -586,18 +595,10 @@ MODULE moduleMesh2DCyl
       !Computes total volume of the cell
       self%volume  = r*detJ*PI8 !4*2*pi
       !Computes volume per node
-      Xi = (/-5.D-1, -5.D-1, 0.D0/)
-      r = self%gatherF(Xi, 4, self%r)
-      self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8
-      Xi = (/ 5.D-1, -5.D-1, 0.D0/)
-      r = self%gatherF(Xi, 4, self%r)
-      self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8
-      Xi = (/ 5.D-1,  5.D-1, 0.D0/)
-      r = self%gatherF(Xi, 4, self%r)
-      self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8
-      Xi = (/-5.D-1,  5.D-1, 0.D0/)
-      r = self%gatherF(Xi, 4, self%r)
-      self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8
+      self%n1%v = self%n1%v + fPsi(1)*self%volume
+      self%n2%v = self%n2%v + fPsi(2)*self%volume
+      self%n3%v = self%n3%v + fPsi(3)*self%volume
+      self%n4%v = self%n4%v + fPsi(4)*self%volume
 
     END SUBROUTINE volumeQuad
 
diff --git a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90 b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
index 2c02358..0df1289 100644
--- a/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
+++ b/src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90
@@ -296,20 +296,6 @@ MODULE moduleMeshInputGmsh2
 
       CLOSE(10)
 
-      ! Adjust node volume at axis
-      SELECT CASE(self%geometry)
-      CASE("Cyl")
-        DO n  = 1, self%numNodes
-          r = self%nodes(n)%obj%getCoordinates()
-          IF (r(2) == 0.D0) THEN
-            self%nodes(n)%obj%v = self%nodes(n)%obj%v * 2.0D0!2.0D0/3.0D0
-
-          END IF
-
-        END DO
-
-      END SELECT
-
       !Call mesh connectivity
       CALL self%connectMesh
 
diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
index eb6af79..309ca5c 100644
--- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
@@ -496,19 +496,19 @@ MODULE moduleMeshInputVTU
 
       END DO
 
-      ! Adjust node volume at axis
-      SELECT CASE(self%geometry)
-      CASE("Cyl")
-        DO n  = 1, numNodes
-          r = self%nodes(n)%obj%getCoordinates()
-          IF (r(2) == 0.D0) THEN
-            self%nodes(n)%obj%v = self%nodes(n)%obj%v * 2.0D0!2.0D0/3.0D0
-
-          END IF
-
-        END DO
-
-      END SELECT
+      ! ! Adjust node volume at axis
+      ! SELECT CASE(self%geometry)
+      ! CASE("Cyl")
+      !   DO n  = 1, numNodes
+      !     r = self%nodes(n)%obj%getCoordinates()
+      !     IF (r(2) == 0.D0) THEN
+      !       self%nodes(n)%obj%v = self%nodes(n)%obj%v * 2.0D0!2.0D0/3.0D0
+      !
+      !     END IF
+      !
+      !   END DO
+      !
+      ! END SELECT
         
       !Call mesh connectivity
       CALL self%connectMesh
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index e96ff2a..af1cd36 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -78,6 +78,8 @@ MODULE moduleMesh
     REAL(8):: normal(1:3)
     !Weight for random injection of particles
     REAL(8):: weight = 1.D0
+    ! Surface of edge
+    REAL(8):: surface = 0.D0
     !Pointer to boundary type
     TYPE(boundaryCont), POINTER:: boundary
     !Array of functions for boundary conditions
diff --git a/src/modules/moduleInject.f90 b/src/modules/moduleInject.f90
index 69a3dd6..1d4e2e7 100644
--- a/src/modules/moduleInject.f90
+++ b/src/modules/moduleInject.f90
@@ -63,6 +63,7 @@ MODULE moduleInject
     INTEGER, ALLOCATABLE:: edges(:) !Array with edges
     REAL(8), ALLOCATABLE:: cumWeight(:) !Array of cummulative probability
     REAL(8):: sumWeight
+    REAL(8):: surface ! Total surface of injection
     TYPE(velDistCont):: v(1:3) !Velocity distribution function in each direction
     CONTAINS
       PROCEDURE, PASS:: init         => initInject
@@ -164,15 +165,12 @@ MODULE moduleInject
 
       END DO
 
-      !Calculates cumulative probability
-      ALLOCATE(self%cumWeight(1:self%nEdges))
-      et = 1
-      self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight
-      DO et = 2, self%nEdges
-        self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1)
+      !Calculates total area
+      self%surface = 0.D0
+      DO et = 1, self%nEdges
+        self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface
 
       END DO
-      self%sumWeight = self%cumWeight(self%nEdges)
 
     END SUBROUTINE initInject
 
@@ -313,12 +311,12 @@ MODULE moduleInject
 
       !$OMP DO
       DO n = nMin, nMax
-        randomX    =  randomWeighted(self%cumWeight, self%sumWeight)
+        randomX    =  random(1, self%nEdges)
         randomEdge => mesh%edges(self%edges(randomX))%obj
         !Random position in edge
         partInj(n)%r = randomEdge%randPos()
         !Assign weight to particle.
-        partInj(n)%weight = self%species%weight
+        partInj(n)%weight = self%species%weight * randomEdge%surface / self%surface
         !Volume associated to the edge:
         IF (ASSOCIATED(randomEdge%e1)) THEN
           partInj(n)%cell = randomEdge%e1%n