Working, but I think I need to redo the injection module to make it more flexible

src/modules/moduleInject.f90

@@ -28,7 +28,7 @@ MODULE moduleInject
   END TYPE injectGeneric
 
   INTEGER:: nInject
-  TYPE(injectGeneric), ALLOCATABLE:: inject(:)
+  TYPE(injectGeneric), ALLOCATABLE, target:: inject(:)
 
   CONTAINS
     !Initialize an injection of particles
@@ -188,14 +188,67 @@ MODULE moduleInject
 
     ! Update the value of injects
     subroutine updateInjects()
+      use moduleMesh, only: meshEdge, meshNode, mesh, qSpecies
       implicit none
       
-      integer:: i
+      integer:: i, e, s, n
+      class(injectGeneric), pointer:: self
+      integer, allocatable:: nodes(:)
+      class(meshEdge), pointer:: edge
+      real(8), allocatable:: density_nodes(:)
+      class(meshNode), pointer:: node
+      real(8):: den_center
+      real(8):: density_incident, density_rest
+      real(8):: alpha
 
       do i = 1, nInject
+        self => inject(i)
         select case(inject(i)%type)
         case ('quasiNeutral')
-          print*, "Calculating quasi-neutrality"
+          do e = 1, self%nEdges
+            edge => self%edges(e)%obj
+
+            density_incident = 0.d0
+            density_rest     = 0.d0
+
+            nodes = edge%getNodes(edge%nNodes)
+            allocate(density_nodes(1:edge%nNodes))
+            do s = 1, nSpecies
+              do n = 1, edge%nNodes
+                node => mesh%nodes(nodes(n))%obj
+
+                density_nodes(n) = node%output(s)%den
+
+              end do
+
+              ! Gather the density at the edge center and multiply by the species charge
+              den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, density_nodes)
+
+              if (s == self%species%n) then
+                density_incident = den_center
+
+              else
+                density_rest = density_rest + den_center
+
+              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
+
+            ! Adjust the weight of particles to match the new current
+            self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-2 * alpha
+
+          end do
 
         end select
       end do