Changed the injection quasiNuetral type to keep a 0 electric field at the edge

src/modules/moduleInject.f90

@@ -53,6 +53,7 @@ MODULE moduleInject
   end type injectconstant
 
   type, extends(injectGeneric):: injectQuasiNeutral
+    real(8):: qSign_incident
 
   end type injectQuasiNeutral
 
@@ -228,6 +229,14 @@ MODULE moduleInject
       !Scale particles for different species steps
       IF (self%nParticles == 0) CALL criticalError("The number of particles for inject is 0.", 'initInject')
 
+      ! Specific parameters per type
+      select type(self)
+      type is(injectQuasiNeutral)
+        ! Get the sign of the charged species
+        self%qSign_incident = sign(1.d0, self%species%qm)
+
+      end select
+
     END SUBROUTINE initInject
 
     ! Update the value of injects
@@ -248,17 +257,15 @@ MODULE moduleInject
 
     ! Updates the flow in the injection to maintain quasineutrality
     subroutine updateQuasiNeutral(self)
-      use moduleMesh, only: meshEdge, meshNode, mesh, qSpecies
+      use moduleMesh, only: meshEdge, meshCell
       implicit none
 
       class(injectGeneric), intent(inout):: self
-      integer:: e, s, n
+      integer:: e
-      integer, allocatable:: nodes(:)
       class(meshEdge), pointer:: edge
-      real(8), allocatable:: density_nodes(:)
+      class(meshCell), pointer:: cell
-      class(meshNode), pointer:: node
+      real(8):: Xi(1:3)
-      real(8):: den_center
+      real(8):: EF_normal
-      real(8):: density_incident, density_rest
       real(8):: alpha
 
       select type(self)
@@ -266,43 +273,22 @@ MODULE moduleInject
         do e = 1, self%nEdges
           edge => self%edges(e)%obj
 
-          density_incident = 0.d0
+          Xi = edge%centerXi()
-          density_rest     = 0.d0
 
-          nodes = edge%getNodes(edge%nNodes)
+          if (associated(edge%e1)) then
-          allocate(density_nodes(1:edge%nNodes))
+            cell => edge%e1
-          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
+            cell => edge%e2
-            alpha = 0.d0
 
           end if
 
+          ! Projection of EF on the edge normal vector
+          EF_normal = dot_product(cell%gatherElectricField(Xi), edge%normal)
+
+          ! Correction for this time step
+          alpha = self%qSign_incident * EF_normal
+
           ! Adjust the weight of particles to match the new current
           self%weightPerEdge(e) = self%weightPerEdge(e) + 1.0d-1 * alpha