Shifting towards constant number of particles per edge

So now each edge has the same number of particles and the weight of each
particle is calculated based on the surface of each edge compared to the
total one.

Only in 2DCyl, still to extend to other geometries.

Not perfect constant density, but the issue might be the node volume.

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/)
+      self%n1%v = self%n1%v + fPsi(1)*self%volume
-      r = self%gatherF(Xi, 4, self%r)
+      self%n2%v = self%n2%v + fPsi(2)*self%volume
-      self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8
+      self%n3%v = self%n3%v + fPsi(3)*self%volume
-      Xi = (/ 5.D-1, -5.D-1, 0.D0/)
+      self%n4%v = self%n4%v + fPsi(4)*self%volume
-      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
 
     END SUBROUTINE volumeQuad
 

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
 

src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90

@@ -496,19 +496,19 @@ MODULE moduleMeshInputVTU
 
       END DO
 
-      ! Adjust node volume at axis
+      ! ! Adjust node volume at axis
-      SELECT CASE(self%geometry)
+      ! SELECT CASE(self%geometry)
-      CASE("Cyl")
+      ! CASE("Cyl")
-        DO n  = 1, numNodes
+      !   DO n  = 1, numNodes
-          r = self%nodes(n)%obj%getCoordinates()
+      !     r = self%nodes(n)%obj%getCoordinates()
-          IF (r(2) == 0.D0) THEN
+      !     IF (r(2) == 0.D0) THEN
-            self%nodes(n)%obj%v = self%nodes(n)%obj%v * 2.0D0!2.0D0/3.0D0
+      !       self%nodes(n)%obj%v = self%nodes(n)%obj%v * 2.0D0!2.0D0/3.0D0
-
+      !
-          END IF
+      !     END IF
-
+      !
-        END DO
+      !   END DO
-
+      !
-      END SELECT
+      ! END SELECT
         
       !Call mesh connectivity
       CALL self%connectMesh

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

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
+      !Calculates total area
-      ALLOCATE(self%cumWeight(1:self%nEdges))
+      self%surface = 0.D0
-      et = 1
+      DO et = 1, self%nEdges
-      self%cumWeight(1) = mesh%edges(self%edges(et))%obj%weight
+        self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface
-      DO et = 2, self%nEdges
-        self%cumWeight(et) = mesh%edges(self%edges(et))%obj%weight + self%cumWeight(et-1)
 
       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