Inject fixed. Now getting close to input

src/modules/init/moduleInput.f90

@@ -1128,12 +1128,35 @@ MODULE moduleInput
       !Builds the K matrix for the Particles mesh
       CALL mesh%constructGlobalK()
 
-      !Assign the procedure to find a volume for meshColl
+      !Assign the procedure to find a cell for meshColl
       IF (doubleMesh) THEN
         findCellColl => findCellCollMesh
+        ! Link edges with cells in meshColl
+        DO e=1, mesh%numEdges
+          nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos())
+          IF (nVolColl > 0) THEN
+            mesh%edges(e)%obj%eColl => meshColl%cells(nVolColl)%obj
+
+          ELSE
+            CALL criticalError("No connection between edge and meshColl", "initInject")
+
+          END IF
+
+        end do
 
       ELSE
         findCellColl => findCellSameMesh
+        ! Link edges with cells in meshColl
+        DO e=1, mesh%numEdges
+          IF (ASSOCIATED(mesh%edges(e)%obj%e1)) THEN
+            mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e1
+
+          ELSE
+            mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e2
+
+          END IF
+
+        end do
 
       END IF
 
@@ -1144,8 +1167,33 @@ MODULE moduleInput
         type is(boundaryQuasiNeutrality)
           ! Loop over all physical surfaces
           do ps = 1, nPhysicalSurfaces
+            ! Loop over all species
             do s = 1, nSpecies
+              ! If the boundary for the species is linked to the one analysing, add the edges
               if (associated(physicalSurfaces(ps)%particles(s), bound)) then
+                bound%edges = [bound%edges, physicalSurfaces(ps)%edges]
+              end if
+
+            end do
+
+          end do
+
+        end select
+
+      end do
+
+      ! EM Boundaries
+      do b = 1, nBoundariesEM
+        select type(bound => boundariesEM(b)%obj)
+        type is(boundaryEMDirichlet)
+          ! Loop over all physical surfaces
+          do ps = 1, nPhysicalSurfaces
+            ! If the boundary for the species is linked to the one analysing, add the edges
+            if (associated(physicalSurfaces(ps)%EM, bound)) then
+              bound%nodes = [bound%nodes, physicalSurfaces(ps)%nodes]
+            end if
+
+          end do
 
         end select
 

src/modules/moduleInject.f90

@@ -2,6 +2,7 @@
 MODULE moduleInject
   USE moduleSpecies
   use velocityDistribution
+  use moduleMesh, only: meshEdgePointer
 
   !Generic injection of particles
   TYPE:: injectGeneric
@@ -14,7 +15,7 @@ MODULE moduleInject
     INTEGER:: nParticles !Number of particles to introduce each time step
     CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
-    INTEGER, ALLOCATABLE:: edges(:) !Array with edges
+    type(meshEdgePointer), allocatable:: edges(:)
     INTEGER, ALLOCATABLE:: particlesPerEdge(:) ! Particles per edge
     REAL(8), ALLOCATABLE:: weightPerEdge(:) ! Weight per edge
     REAL(8):: surface ! Total surface of injection
@@ -43,6 +44,7 @@ MODULE moduleInject
       INTEGER, INTENT(in):: i
       REAL(8), INTENT(in):: v, n(1:3), temperature(1:3)
       INTEGER, INTENT(in):: sp, ps, particlesPerEdge
+      integer:: ps_index
       REAL(8):: tauInject
       REAL(8), INTENT(in):: flow
       CHARACTER(:), ALLOCATABLE, INTENT(in):: units
@@ -55,50 +57,21 @@ MODULE moduleInject
       self%vMod        =  v / v_ref
       self%n           =  n / NORM2(n)
       self%temperature =  temperature / T_ref
-      !Gets the edge elements from which particles are injected
-      DO e = 1, mesh%numEdges
-        phSurface(e) = mesh%edges(e)%obj%physicalSurface
 
-      END DO
+      !Get array index of physical surface for inject
-      self%nEdges = COUNT(phSurface == physicalSurface)
+      ps_index = physicalSurface_to_index(ps)
-      ALLOCATE(self%edges(1:self%nEdges))
-      ALLOCATE(self%particlesPerEdge(1:self%nEdges))
-      ALLOCATE(self%weightPerEdge(1:self%nEdges))
-      et = 0
-      DO e=1, mesh%numEdges
-        IF (mesh%edges(e)%obj%physicalSurface == physicalSurface) THEN
-          et = et + 1
-          self%edges(et) = mesh%edges(e)%obj%n
-          !Assign connectivity between injection edge and meshColl volume
-          IF (doubleMesh) THEN
-            nVolColl = findCellBrute(meshColl, mesh%edges(e)%obj%randPos())
-            IF (nVolColl > 0) THEN
-              mesh%edges(e)%obj%eColl => meshColl%cells(nVolColl)%obj
 
-            ELSE
+      ! Copy array of edges
-              CALL criticalError("No connection between edge and meshColl", "initInject")
+      self%edges = physicalSurfaces(ps_index)%edges
+      self%nEdges = size(self%edges)
 
-            END IF
+      allocate(self%particlesPerEdge(1:self%nEdges))
-
+      allocate(self%weightPerEdge(1:self%nEdges))
-          ELSE
-            IF (ASSOCIATED(mesh%edges(e)%obj%e1)) THEN
-              mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e1
-
-            ELSE
-              mesh%edges(e)%obj%eColl => mesh%edges(e)%obj%e2
-
-            END IF
-
-          END IF
-
-        END IF
-
-      END DO
 
       !Calculates total area
       self%surface = 0.D0
       DO et = 1, self%nEdges
-        self%surface = self%surface + mesh%edges(self%edges(et))%obj%surface
+        self%surface = self%surface + self%edges(et)%obj%surface
 
       END DO
 
@@ -149,7 +122,7 @@ MODULE moduleInject
         ! Particles per edge defined by the user
         self%particlesPerEdge = particlesPerEdge
         DO et = 1, self%nEdges
-          self%weightPerEdge(et) = fluxPerStep*mesh%edges(self%edges(et))%obj%surface / REAL(particlesPerEdge)
+          self%weightPerEdge(et) = fluxPerStep * self%edges(et)%obj%surface / real(particlesPerEdge)
 
         END DO
 
@@ -159,7 +132,7 @@ MODULE moduleInject
         ! No particles assigned per edge, use the species weight
         self%weightPerEdge = self%species%weight
         DO et = 1, self%nEdges
-          self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*mesh%edges(self%edges(et))%obj%surface / self%species%weight))
+          self%particlesPerEdge(et) = max(1,FLOOR(fluxPerStep*self%edges(et)%obj%surface / self%species%weight))
         END DO
 
         self%nParticles = SUM(self%particlesPerEdge)
@@ -266,7 +239,7 @@ MODULE moduleInject
       !$OMP DO
       DO e = 1, self%nEdges
         ! Select edge for injection
-        randomEdge => mesh%edges(self%edges(e))%obj
+        randomEdge => self%edges(e)%obj
         ! Inject particles in edge
         DO i = 1, self%particlesPerEdge(e)
           ! Index in the global partInj array