Cylindrical injection working better

Seems things are a bit better. Still, more cases are needed and still
not perfectly uniform...

src/modules/common/moduleRandom.f90

@@ -75,10 +75,10 @@ MODULE moduleRandom
       REAL(8):: rnd0b
       INTEGER:: rnd, i
 
-      rnd0b = random(0.D0, sumWeight)
+      rnd0b = random()
       i = 1
       DO 
-        IF (rnd0b <= cumWeight(i)) THEN
+        IF (rnd0b <= cumWeight(i)/sumWeight) THEN
           rnd = i
           EXIT
 

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

@@ -163,7 +163,7 @@ MODULE moduleMesh2DCyl
       r2 = self%n2%getCoordinates()
       self%z  = (/r1(1), r2(1)/)
       self%r  = (/r1(2), r2(2)/)
-      self%weight = r2(2)**2 - r1(2)**2
+      self%weight = DABS(self%r(2)**2 - self%r(1)**2)
       !Normal vector
       self%normal = (/ -(self%r(2)-self%r(1)), &
                          self%z(2)-self%z(1) , &
@@ -223,21 +223,12 @@ MODULE moduleMesh2DCyl
       CLASS(meshEdge2DCyl), INTENT(in):: self
       REAL(8):: rnd
       REAL(8):: r(1:3)
-      REAL(8):: dr, dz
+      REAL(8):: p1(1:2), p2(1:2)
 
       rnd = random()
-      dr  = self%r(2) - self%r(1)
-      dz  = self%z(2) - self%z(1)
-      IF (dr /= 0.D0) THEN
-        r(2) = dr * DSQRT(rnd) + self%r(1)
-        r(1) = dz * (r(2) - self%r(1))/dr + self%z(1)
-
-      ELSE
-        r(2) = self%r(1)
-        r(1) = dz * rnd + self%z(1)
-
-      END IF
-
+      p1 = (/self%z(1), self%r(1) /)
+      p2 = (/self%z(2), self%r(2) /)
+      r(1:2) = (1.D0 - rnd)*p1 + rnd*p2
       r(3)   =  0.D0
 
     END FUNCTION randPosEdge
@@ -572,7 +563,7 @@ MODULE moduleMesh2DCyl
 
     !Compute element volume
     PURE SUBROUTINE volumeQuad(self)
-      USE moduleConstParam, ONLY: PI8
+      USE moduleConstParam, ONLY: PI, PI8
       IMPLICIT NONE
 
       CLASS(meshCell2DCylQuad), INTENT(inout):: self
@@ -596,15 +587,31 @@ MODULE moduleMesh2DCyl
       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
+      IF (self%r(1) == 0.D0) THEN
+        self%n1%v = self%n1%v * PI / 2.D0
+
+      END IF
       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
+      IF (self%r(2) == 0.D0) THEN
+        self%n2%v = self%n2%v * PI / 2.D0
+
+      END IF
       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
+      IF (self%r(3) == 0.D0) THEN
+        self%n3%v = self%n3%v * PI / 2.D0
+
+      END IF
       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
+      IF (self%r(4) == 0.D0) THEN
+        self%n4%v = self%n4%v * PI / 2.D0
+
+      END IF
 
     END SUBROUTINE volumeQuad