Finally, using the right center for each element

src/modules/mesh/1DCart/moduleMesh1DCart.f90

@@ -29,6 +29,7 @@ MODULE moduleMesh1DCart
       PROCEDURE, PASS::   intersection => intersection1DCart
       PROCEDURE, PASS::   randPos      => randPosEdge
       procedure, pass::   center       => centerEdgePoint
+      procedure, nopass:: centerXi     => centerXiPoint
       procedure, nopass:: fPsi         => fPsiPoint
 
   END TYPE meshEdge1DCart

src/modules/mesh/1DRad/moduleMesh1DRad.f90

@@ -29,6 +29,7 @@ MODULE moduleMesh1DRad
       PROCEDURE, PASS::   intersection => intersection1DRad
       PROCEDURE, PASS::   randPos      => randPos1DRad
       procedure, pass::   center       => centerEdgePoint
+      procedure, nopass:: centerXi     => centerXiPoint
       procedure, nopass:: fPsi         => fPsiPoint
 
   END TYPE meshEdge1DRad

src/modules/mesh/2DCart/moduleMesh2DCart.f90

@@ -29,6 +29,7 @@ MODULE moduleMesh2DCart
       PROCEDURE, PASS::   intersection => intersection2DCartEdge
       PROCEDURE, PASS::   randPos      => randPosEdge
       procedure, pass::   center       => centerEdgeSegm
+      procedure, nopass:: centerXi     => centerXiSegm
       procedure, nopass:: fPsi         => fPsiSegm
 
   END TYPE meshEdge2DCart

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

@@ -29,6 +29,7 @@ MODULE moduleMesh2DCyl
       PROCEDURE, PASS::   intersection => intersection2DCylEdge
       PROCEDURE, PASS::   randPos      => randPosEdge
       procedure, pass::   center       => centerEdgeSegm
+      procedure, nopass:: centerXi     => centerXiSegm
       procedure, nopass:: fPsi         => fPsiSegm
 
   END TYPE meshEdge2DCyl

src/modules/mesh/3DCart/moduleMesh3DCart.f90

@@ -31,6 +31,7 @@ MODULE moduleMesh3DCart
       PROCEDURE, PASS:: randPos       => randPosEdgeTria
       procedure, pass:: center        => centerEdgeTria
       !PARTICULAR PROCEDURES
+      procedure, nopass:: centerXi => centerXiTria
       PROCEDURE, NOPASS:: fPsi     => fPsiTria
 
   END TYPE meshEdge3DCartTria

src/modules/mesh/moduleMesh.f90

@@ -119,6 +119,7 @@ MODULE moduleMesh
       PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection
       PROCEDURE(randPosEdge_interface),      DEFERRED, PASS:: randPos
       procedure(centerEdge_interface),       deferred, pass:: center
+      procedure(centerXiEdge_interface),     deferred, nopass:: centerXi
       PROCEDURE(fPsi_interface),             DEFERRED, NOPASS:: fPsi
       !Gather value and spatial derivative on the nodes at position Xi
       PROCEDURE, PASS, PRIVATE:: gatherF_edge_scalar
@@ -185,6 +186,12 @@ MODULE moduleMesh
 
     END FUNCTION centerEdge_interface
 
+    ! Returns the center point of an edge in natural coordinates
+    FUNCTION centerXiEdge_interface() RESULT(Xi)
+      REAL(8):: Xi(1:3)
+
+    END FUNCTION centerXiEdge_interface
+
   END INTERFACE
 
   !Containers for edges in the mesh

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -346,7 +346,6 @@ submodule(moduleMesh) boundaryParticle
 
       class(boundaryParticleGeneric), intent(inout):: self
       integer, intent(in):: fileID
-      integer:: e
 
       write(fileID, '(A)') self%name
       select type(self)
@@ -432,18 +431,26 @@ submodule(moduleMesh) boundaryParticle
             end do
 
             if (s == self%s_incident) then
-              density_incident = edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes)
+              density_incident = edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
 
             else
-              density_rest = density_rest + edge%gatherF((/0.d0,0.d0,0.d0/), edge%nNodes, density_nodes)
+              density_rest = density_rest + edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
 
             end if
 
           end do
 
           ! Correction for this time step
+          if (density_rest > 1.0d-10) then
+            ! If there is a rest population, correct
             alpha = 1.d0 - density_incident/density_rest
 
+          else
+            ! If not, alpha is assumed unchaged
+            alpha = 0.d0
+
+          end if
+
           ! Apply correction with a factor of 0.1 to avoid fast changes
           self%alpha(edge%n) = self%alpha(edge%n) + 1.0d-2 * alpha
 
@@ -472,7 +479,7 @@ submodule(moduleMesh) boundaryParticle
       type is(boundaryQuasiNeutrality)
         write(fileID, '(A,",",A)') '"Edge id"', '"alpha"'
         do e = 1, size(self%edges)
-          write(fileID, '('//fmtInt//',",",'//fmtReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n)
+          write(fileID, '('//fmtColInt//','//fmtReal//')') self%edges(e)%obj%n, self%alpha(self%edges(e)%obj%n)
 
         end do
 

src/modules/mesh/moduleMeshCommon.f90

@@ -15,14 +15,69 @@ module moduleMeshCommon
   real(8), parameter:: wQuad(1:3)   = (/        5.D0/9.D0,   8.D0/9.D0,       5.D0/9.D0  /)
 
   ! Center point in natural coordinates
+  ! Point
+  real(8), parameter:: cenPoint(1:3) = (/ 0.0d0,       0.0d0,       0.0d0 /)
   ! Segment
-  real(8), parameter:: cenSeg(1:3)  = (/ 0.5d0,       0.0d0,       0.0d0 /)
+  real(8), parameter:: cenSeg(1:3)   = (/ 0.0d0,       0.0d0,       0.0d0 /)
   ! Tria
   real(8), parameter:: cenTria(1:3)  = (/ 1.0d0/3.0d0, 1.0d0/3.0d0, 0.0d0 /)
   ! Quad
   real(8), parameter:: cenQuad(1:3)  = (/ 0.0d0,       0.0d0,       0.0d0 /)
+  ! Tetra
+  real(8), parameter:: cenTetra(1:3) = (/ 1.0d0/3.0d0, 1.0d0/3.0d0, 1.d0/3.d0 /)
 
   contains
+    ! RETURN CENTER IN NATURAL COORDINATES
+    ! Point
+    pure function centerXiPoint() result(Xi)
+      implicit none
+
+      real(8):: Xi(1:3)
+
+      Xi = cenPoint
+
+    end function centerXiPoint
+
+    ! Segment
+    pure function centerXiSegm() result(Xi)
+      implicit none
+
+      real(8):: Xi(1:3)
+
+      Xi = cenSeg
+
+    end function centerXiSegm
+
+    ! Tria
+    pure function centerXiTria() result(Xi)
+      implicit none
+
+      real(8):: Xi(1:3)
+
+      Xi = cenTria
+
+    end function centerXiTria
+
+    ! Quad
+    pure function centerXiQuad() result(Xi)
+      implicit none
+
+      real(8):: Xi(1:3)
+
+      Xi = cenQuad
+
+    end function centerXiQuad
+
+    ! Quad
+    pure function centerXiTetra() result(Xi)
+      implicit none
+
+      real(8):: Xi(1:3)
+
+      Xi = cenTetra
+
+    end function centerXiTetra
+
     ! ELEMENT FUNCTIONS
     ! Point
     pure function fPsiPoint(Xi, nNodes) RESULT(fPsi)