Add new elements to point containers is now more general

src/modules/init/moduleInput.f90

@@ -1137,6 +1137,20 @@ MODULE moduleInput
 
       END IF
 
+      ! If needed, add nodes and edges to boundary models
+      ! Particle boundaries
+      do b = 1, nBoundariesParticle
+        select type(bound => boundariesParticle(b)%obj)
+        type is(boundaryQuasiNeutrality)
+        ! Loop over all physical surfaces
+        do ps = 1, nPhysicalSurfaces
+          do s = 1, nSpecies
+            if (associated(physicalSurfaces(ps)%particles(s), bound)) then
+
+        end select
+
+      end do
+
     END SUBROUTINE readGeometry
 
     SUBROUTINE readProbes(config)

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

@@ -127,8 +127,8 @@ MODULE moduleMesh1DCart
       ps_index = physicalSurface_to_index(ps)
 
       ! Add elements to physical surface
-      call physicalSurfaces(ps_index)%addEdge(self%n)
-      call physicalSurfaces(ps_index)%addNode(self%n1%n)
+      call meshEdgePointer_add(physicalSurfaces(ps_index)%edges, self%n)
+      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n1%n)
 
     end subroutine initEdge1DCart
 

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

@@ -127,8 +127,8 @@ MODULE moduleMesh1DRad
       ps_index = physicalSurface_to_index(ps)
 
       ! Add elements to physical surface
-      call physicalSurfaces(ps_index)%addEdge(self%n)
-      call physicalSurfaces(ps_index)%addNode(self%n1%n)
+      call meshEdgePointer_add(physicalSurfaces(ps_index)%edges, self%n)
+      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n1%n)
 
     end subroutine initEdge1DRad
 

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

@@ -166,9 +166,9 @@ MODULE moduleMesh2DCart
       ps_index = physicalSurface_to_index(ps)
 
       ! Add elements to physical surface
-      call physicalSurfaces(ps_index)%addEdge(self%n)
-      call physicalSurfaces(ps_index)%addNode(self%n1%n)
-      call physicalSurfaces(ps_index)%addNode(self%n2%n)
+      call meshEdgePointer_add(physicalSurfaces(ps_index)%edges, self%n)
+      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n1%n)
+      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n2%n)
 
     END SUBROUTINE initEdge2DCart
 

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

@@ -175,9 +175,9 @@ MODULE moduleMesh2DCyl
       ps_index = physicalSurface_to_index(ps)
 
       ! Add elements to physical surface
-      call physicalSurfaces(ps_index)%addEdge(self%n)
-      call physicalSurfaces(ps_index)%addNode(self%n1%n)
-      call physicalSurfaces(ps_index)%addNode(self%n2%n)
+      call meshEdgePointer_add(physicalSurfaces(ps_index)%edges, self%n)
+      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n1%n)
+      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n2%n)
 
     END SUBROUTINE initEdge2DCyl
 

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

@@ -148,10 +148,10 @@ MODULE moduleMesh3DCart
       ps_index = physicalSurface_to_index(ps)
 
       ! Add elements to physical surface
-      call physicalSurfaces(ps_index)%addEdge(self%n)
-      call physicalSurfaces(ps_index)%addNode(self%n1%n)
-      call physicalSurfaces(ps_index)%addNode(self%n2%n)
-      call physicalSurfaces(ps_index)%addNode(self%n3%n)
+      call meshEdgePointer_add(physicalSurfaces(ps_index)%edges, self%n)
+      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n1%n)
+      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n2%n)
+      call meshNodePointer_add(physicalSurfaces(ps_index)%nodes, self%n3%n)
 
     END SUBROUTINE initEdge3DCartTria
 

src/modules/mesh/moduleMesh.f90

@@ -64,7 +64,6 @@ MODULE moduleMesh
   !Containers for nodes in the mesh
   TYPE:: meshNodeCont
     CLASS(meshNode), ALLOCATABLE:: obj
-    CONTAINS
 
   END TYPE meshNodeCont
 
@@ -79,6 +78,12 @@ MODULE moduleMesh
   END TYPE meshNodePointer
 
   interface
+    module subroutine meshNodePointer_add(self, node)
+      class(meshNodePointer), allocatable, intent(inout):: self(:)
+      integer, intent(in):: node
+
+    end subroutine meshNodePointer_add
+
     module function meshNodePointer_equal_type_type(self, other) result(isEqual)
       class(meshNodePointer), intent(in):: self, other
       logical:: isEqual
@@ -192,6 +197,12 @@ MODULE moduleMesh
   end type meshEdgePointer
 
   interface
+    module subroutine meshEdgePointer_add(self, edge)
+      class(meshEdgePointer), allocatable, intent(inout):: self(:)
+      integer, intent(in):: edge
+
+    end subroutine meshEdgePointer_add
+
     module function meshEdgePointer_equal_type_type(self, other) result(isEqual)
       class(meshEdgePointer), intent(in):: self, other
       logical:: isEqual
@@ -837,7 +848,7 @@ MODULE moduleMesh
   end type boundaryParticlePointer
 
   !Number of boundaries
-  INTEGER:: nBoundaries = 0
+  INTEGER:: nBoundariesParticle = 0
   !Array for boundaries
   TYPE(boundaryParticleCont), ALLOCATABLE, TARGET:: boundariesParticle(:)
 
@@ -853,7 +864,6 @@ MODULE moduleMesh
     contains
       procedure, pass::                              init   => initBoundaryEM
       procedure(applyEM_interface), deferred, pass:: apply
-      procedure, pass::                              addNodes
 
   end type boundaryEMGeneric
 
@@ -868,12 +878,6 @@ MODULE moduleMesh
 
     end subroutine initBoundaryEM
 
-    module subroutine addNodes(self, nodes)
-      class(boundaryEMGeneric), intent(inout):: self
-      integer, intent(in):: nodes(:)
-
-    end subroutine addNodes
-
     module function boundaryEMName_to_Index(boundaryName) result(bp)
       character(:), allocatable:: boundaryName
       integer:: bp
@@ -957,27 +961,9 @@ MODULE moduleMesh
     class(meshEdgePointer), allocatable:: edges(:)
     class(boundaryParticlePointer), allocatable:: particles(:)
     class(boundaryEMGeneric), pointer:: EM => null()
-    contains
-      procedure, pass:: addNode
-      procedure, pass:: addEdge
 
   end type
 
-  interface
-    module subroutine addNode(self, node)
-      class(physicalSurface), intent(inout):: self
-      integer, intent(in):: node
-
-    end subroutine addNode
-
-    module subroutine addEdge(self, edge)
-      class(physicalSurface), intent(inout):: self
-      integer, intent(in):: edge
-
-    end subroutine addEdge
-
-  end interface
-
   integer:: nPhysicalSurfaces
   type(physicalSurface), allocatable:: physicalSurfaces(:)
 

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -9,7 +9,7 @@ submodule(moduleMesh) boundaryEM
       integer:: b
 
       bp = 0
-      do b = 1, nBoundaries
+      do b = 1, nBoundariesEM
         if (boundaryName == boundariesEM(b)%obj%name) then
           bp = boundariesEM(b)%obj%n
 
@@ -37,6 +37,7 @@ submodule(moduleMesh) boundaryEM
       logical:: found
 
       self%n = b
+      allocate(self%nodes(0))
       call config%get(object // '.name', self%name, found)
       if (.not. found) then
         call criticalError('Required parameter "name" for EM boundary condition not found', &

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -10,7 +10,7 @@ submodule(moduleMesh) boundaryParticle
       integer:: b
 
       bp = 0
-      do b = 1, nBoundaries
+      do b = 1, nBoundariesParticle
         if (boundaryName == boundariesParticle(b)%obj%name) then
           bp = boundariesParticle(b)%obj%n
 

src/modules/mesh/moduleMesh@elements.f90

@@ -18,6 +18,37 @@ submodule(moduleMesh) elements
 
     END SUBROUTINE resetOutput
 
+    subroutine meshNodePointerAdd(self, node)
+      implicit none
+
+      class(meshNodePointer), allocatable, intent(inout):: self(:)
+      integer, intent(in):: node
+      integer:: n
+      logical:: inArray
+      type(meshNodePointer):: nodeToAdd
+
+
+      nodeToAdd%obj => mesh%nodes(node)%obj
+      inArray = .false.
+      ! I cannot use the procedure 'any' for this
+      do n = 1 , size(self)
+        IF (self(n) == nodeToAdd) THEN
+          ! Node already in array
+          inArray = .true.
+          exit
+
+        end if
+
+      end do
+
+      if (.not. inArray) then
+        ! If not, add element to array of nodes
+        self = [self, nodeToAdd]
+
+      end if
+
+    end subroutine meshNodePointerAdd
+
     function meshNodePointer_equal_type_type(self, other) result(isEqual)
       implicit none
 
@@ -49,6 +80,37 @@ submodule(moduleMesh) elements
 
     end function meshEdgePointer_equal_type_type
 
+    subroutine meshEdgePointer_add(self, edge)
+      implicit none
+
+      class(meshEdgePointer), allocatable, intent(inout):: self(:)
+      integer, intent(in):: edge
+      integer:: n
+      logical:: inArray
+      type(meshEdgePointer):: edgeToAdd
+
+
+      edgeToAdd%obj => mesh%edges(edge)%obj
+      inArray = .false.
+      ! I cannot use the procedure 'any' for this
+      do n = 1 , size(self)
+        IF (self(n) == edgeToAdd) THEN
+          ! Node already in array
+          inArray = .true.
+          exit
+
+        end if
+
+      end do
+
+      if (.not. inArray) then
+        ! If not, add element to array of nodes
+        self = [self, edgeToAdd]
+
+      end if
+
+    end subroutine meshEdgePointer_add
+
     function meshEdgePointer_equal_type_int(self, other) result(isEqual)
       implicit none
 

src/modules/mesh/moduleMesh@surfaces.f90

@@ -19,66 +19,4 @@ submodule(moduleMesh) surfaces
 
     end function physicalSurface_to_index
 
-    subroutine addNode(self, node)
-      implicit none
-
-      class(physicalSurface), intent(inout):: self
-      integer, intent(in):: node
-      integer:: nn
-      logical:: inArray
-      type(meshNodePointer):: nodeToAdd
-
-
-      nodeToAdd%obj => mesh%nodes(node)%obj
-      inArray = .false.
-      ! I cannot use the procedure 'any' for this
-      do nn = 1 , size(self%nodes)
-        IF (self%nodes(nn) == nodeToAdd) THEN
-          ! Node already in array
-          inArray = .true.
-          exit
-
-        end if
-
-      end do
-
-      if (.not. inArray) then
-        ! If not, add element to array of nodes
-        self%nodes = [self%nodes, nodeToAdd]
-
-      end if
-
-    end subroutine addNode
-
-    subroutine addEdge(self, edge)
-      implicit none
-
-      class(physicalSurface), intent(inout):: self
-      integer, intent(in):: edge
-      integer:: nn
-      logical:: inArray
-      type(meshEdgePointer):: edgeToAdd
-
-
-      edgeToAdd %obj => mesh%edges(edge)%obj
-      inArray = .false.
-      ! I cannot use the procedure 'any' for this
-      do nn = 1 , size(self%edges)
-        IF (self%edges(nn) == edgeToAdd) THEN
-          ! Node already in array
-          inArray = .true.
-          exit
-
-        end if
-
-      end do
-
-      if (.not. inArray) then
-        ! If not, add element to array of nodes
-        self%edges = [self%edges, edgeToAdd]
-
-      end if
-
-    end subroutine addEdge
-
 end submodule surfaces

src/modules/moduleInject.f90

@@ -30,7 +30,7 @@ MODULE moduleInject
 
   CONTAINS
     !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, physicalSurface, particlesPerEdge)
+    SUBROUTINE initInject(self, i, v, n, temperature, flow, units, sp, ps, particlesPerEdge)
       USE moduleMesh
       USE moduleRefParam
       USE moduleConstParam
@@ -42,7 +42,7 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(inout):: self
       INTEGER, INTENT(in):: i
       REAL(8), INTENT(in):: v, n(1:3), temperature(1:3)
-      INTEGER, INTENT(in):: sp, physicalSurface, particlesPerEdge
+      INTEGER, INTENT(in):: sp, ps, particlesPerEdge
       REAL(8):: tauInject
       REAL(8), INTENT(in):: flow
       CHARACTER(:), ALLOCATABLE, INTENT(in):: units
@@ -312,7 +312,7 @@ MODULE moduleInject
             !If while injecting a no-drift distribution the velocity is negative, reflect
             if ((self%vMod == 0.D0) .and. &
                 (dot_product(direction, partInj(n)%v) < 0.D0)) then
-              call reflection(randomEdge, partInj(n))
+              call genericReflection(randomEdge, partInj(n))
 
             end if