Working towards compilation with submodules

src/makefile

@@ -1,7 +1,7 @@
 OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/moduleCompTime.o \
 			  	$(OBJDIR)/moduleSpecies.o $(OBJDIR)/moduleInject.o $(OBJDIR)/moduleInput.o \
 			  	$(OBJDIR)/moduleErrors.o $(OBJDIR)/moduleList.o $(OBJDIR)/moduleOutput.o \
-			  	$(OBJDIR)/moduleBoundary.o $(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
+			  	$(OBJDIR)/moduleCaseParam.o $(OBJDIR)/moduleRefParam.o \
 			  	$(OBJDIR)/moduleCollisions.o $(OBJDIR)/moduleTable.o $(OBJDIR)/moduleParallel.o \
 			  	$(OBJDIR)/moduleEM.o $(OBJDIR)/moduleRandom.o $(OBJDIR)/moduleMath.o \
 					$(OBJDIR)/moduleProbe.o $(OBJDIR)/moduleAverage.o $(OBJDIR)/moduleCoulomb.o \
@@ -9,7 +9,7 @@ OBJECTS = $(OBJDIR)/moduleMesh.o $(OBJDIR)/moduleMeshBoundary.o $(OBJDIR)/module
 			  	$(OBJDIR)/moduleMeshInputVTU.o $(OBJDIR)/moduleMeshOutputVTU.o \
 			  	$(OBJDIR)/moduleMeshInputGmsh2.o $(OBJDIR)/moduleMeshOutputGmsh2.o \
 			  	$(OBJDIR)/moduleMeshInput0D.o $(OBJDIR)/moduleMeshOutput0D.o \
-				$(OBJDIR)/moduleMeshInputText.o $(OBJDIR)/moduleMeshOutputText.o \
+					$(OBJDIR)/moduleMeshInputText.o $(OBJDIR)/moduleMeshOutputText.o \
 			  	$(OBJDIR)/moduleMesh3DCart.o \
 			  	$(OBJDIR)/moduleMesh2DCyl.o  \
 			  	$(OBJDIR)/moduleMesh2DCart.o \

src/modules/init/moduleInput.f90

@@ -8,7 +8,6 @@ MODULE moduleInput
     SUBROUTINE readConfig(inputFile)
       USE json_module
       USE moduleErrors
-      USE moduleBoundary
       USE moduleOutput
       USE moduleMesh
       IMPLICIT NONE
@@ -793,7 +792,7 @@ MODULE moduleInput
 
     !Reads boundary conditions for the mesh
     SUBROUTINE readBoundary(config)
-      USE moduleBoundary
+      use moduleMesh
       USE moduleErrors
       USE moduleSpecies
       USE moduleRefParam
@@ -817,22 +816,22 @@ MODULE moduleInput
       INTEGER:: nTypes
 
       CALL config%info('boundary', found, n_children = nBoundary)
-      ALLOCATE(boundary(1:nBoundary))
+      ALLOCATE(boundaries(1:nBoundary))
       DO i = 1, nBoundary
         WRITE(iString, '(i2)') i
         object = 'boundary(' // TRIM(iString) // ')'
 
-        boundary(i)%n = i
+        boundaries(i)%n = i
-        CALL config%get(object // '.name',            boundary(i)%name,            found)
+        CALL config%get(object // '.name',            boundaries(i)%name,            found)
-        CALL config%get(object // '.physicalSurface', boundary(i)%physicalSurface, found)
+        CALL config%get(object // '.physicalSurface', boundaries(i)%physicalSurface, found)
         CALL config%info(object // '.bTypes', found, n_children = nTypes)
         IF (nTypes /= nSpecies) CALL criticalError('Not enough boundary types defined in ' // object, 'readBoundary')
-        ALLOCATE(boundary(i)%bTypes(1:nSpecies))
+        ALLOCATE(boundaries(i)%bTypes(1:nSpecies))
         DO s = 1, nSpecies
           WRITE(sString,'(i2)') s
           object = 'boundary(' // TRIM(iString) // ').bTypes(' // TRIM(sString) // ')'
           CALL config%get(object // '.type', bType, found)
-          associate(bound => boundary(i)%bTypes(s)%obj)
+          associate(bound => boundaries(i)%bTypes(s)%obj)
             SELECT CASE(bType)
             CASE('reflection')
               ALLOCATE(boundaryReflection:: bound)

src/modules/makefile

@@ -1,5 +1,5 @@
 OBJS = common.o output.o mesh.o solver.o init.o \
-			 moduleBoundary.o moduleCollisions.o moduleInject.o \
+			 moduleCollisions.o moduleInject.o \
 			 moduleList.o moduleProbe.o moduleCoulomb.o \
 			 moduleSpecies.o
 
@@ -11,7 +11,7 @@ common.o:
 output.o: moduleSpecies.o common.o
 	$(MAKE) -C output all
 
-mesh.o: moduleCollisions.o moduleCoulomb.o moduleBoundary.o output.o common.o 
+mesh.o: moduleList.o moduleSpecies.o moduleCollisions.o moduleCoulomb.o output.o common.o 
 	$(MAKE) -C mesh all
 
 solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o
@@ -20,9 +20,6 @@ solver.o: moduleSpecies.o moduleProbe.o common.o output.o mesh.o
 init.o: common.o solver.o moduleInject.o
 	$(MAKE) -C init all
 
-moduleBoundary.o: common.o moduleBoundary.f90
-	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
-
 moduleCollisions.o: moduleList.o moduleSpecies.o common.o moduleCollisions.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 

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

@@ -4,7 +4,6 @@
 !              z == unused
 MODULE moduleMesh1DCart
   USE moduleMesh
-  USE moduleMeshBoundary
   IMPLICIT NONE
 
   REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0,       DSQRT(3.D0/5.D0) /)
@@ -103,7 +102,6 @@ MODULE moduleMesh1DCart
     !Init edge element
     SUBROUTINE initEdge1DCart(self, n, p, bt, physicalSurface)
       USE moduleSpecies
-      USE moduleBoundary
       USE moduleErrors
       IMPLICIT NONE
 
@@ -128,7 +126,7 @@ MODULE moduleMesh1DCart
       self%normal = (/ 1.D0, 0.D0, 0.D0 /)
 
       !Boundary index
-      self%boundary => boundary(bt)
+      self%boundary => boundaries(bt)
       ALLOCATE(self%fboundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies

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

@@ -4,7 +4,6 @@
 !              z == unused
 MODULE moduleMesh1DRad
   USE moduleMesh
-  USE moduleMeshBoundary
   IMPLICIT NONE
 
   REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0,       DSQRT(3.D0/5.D0) /)
@@ -103,7 +102,6 @@ MODULE moduleMesh1DRad
     !Init edge element
     SUBROUTINE initEdge1DRad(self, n, p, bt, physicalSurface)
       USE moduleSpecies
-      USE moduleBoundary
       USE moduleErrors
       IMPLICIT NONE
 
@@ -128,7 +126,7 @@ MODULE moduleMesh1DRad
       self%normal = (/ 1.D0, 0.D0, 0.D0 /)
 
       !Boundary index
-      self%boundary => boundary(bt)
+      self%boundary => boundaries(bt)
       ALLOCATE(self%fboundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies

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

@@ -4,7 +4,6 @@
 !                  z == unused
 MODULE moduleMesh2DCart
   USE moduleMesh
-  USE moduleMeshBoundary
   IMPLICIT NONE
 
   !Values for Gauss integral
@@ -143,7 +142,6 @@ MODULE moduleMesh2DCart
     !Init edge element
     SUBROUTINE initEdge2DCart(self, n, p, bt, physicalSurface)
       USE moduleSpecies
-      USE moduleBoundary
       USE moduleErrors
       IMPLICIT NONE
 
@@ -172,7 +170,7 @@ MODULE moduleMesh2DCart
       self%normal = self%normal/NORM2(self%normal)
 
       !Boundary index
-      self%boundary => boundary(bt)
+      self%boundary => boundaries(bt)
       ALLOCATE(self%fboundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies

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

@@ -4,7 +4,6 @@
 !               z == theta (unused)
 MODULE moduleMesh2DCyl
   USE moduleMesh
-  USE moduleMeshBoundary
   IMPLICIT NONE
 
   !Values for Gauss integral
@@ -143,7 +142,6 @@ MODULE moduleMesh2DCyl
     !Init edge element
     SUBROUTINE initEdge2DCyl(self, n, p, bt, physicalSurface)
       USE moduleSpecies
-      USE moduleBoundary
       USE moduleErrors
       USE moduleConstParam, ONLY: PI
       IMPLICIT NONE
@@ -181,7 +179,7 @@ MODULE moduleMesh2DCyl
       self%normal = self%normal/NORM2(self%normal)
 
       !Boundary index
-      self%boundary => boundary(bt)
+      self%boundary => boundaries(bt)
       ALLOCATE(self%fboundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies

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

@@ -4,7 +4,6 @@
 !                  z == z
 MODULE moduleMesh3DCart
   USE moduleMesh
-  USE moduleMeshBoundary
   IMPLICIT NONE
 
   TYPE, PUBLIC, EXTENDS(meshNode):: meshNode3DCart
@@ -106,7 +105,6 @@ MODULE moduleMesh3DCart
     !Init surface element
     SUBROUTINE initEdge3DCartTria(self, n, p, bt, physicalSurface)
       USE moduleSpecies
-      USE moduleBoundary
       USE moduleErrors
       USE moduleMath
       USE moduleRefParam, ONLY: L_ref
@@ -146,7 +144,7 @@ MODULE moduleMesh3DCart
       self%surface = 1.D0/L_ref**2 !TODO: FIX THIS WHEN MOVING TO 3D
 
       !Boundary index
-      self%boundary => boundary(bt)
+      self%boundary => boundaries(bt)
       ALLOCATE(self%fBoundary(1:nSpecies))
       !Assign functions to boundary
       DO s = 1, nSpecies

src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90

@@ -30,7 +30,6 @@ MODULE moduleMeshInputGmsh2
       USE moduleMesh2DCart
       USE moduleMesh1DRad
       USE moduleMesh1DCart
-      USE moduleBoundary
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(inout):: self

src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90

@@ -161,7 +161,6 @@ MODULE moduleMeshInputVTU
       USE moduleMesh2DCart
       USE moduleMesh1DRad
       USE moduleMesh1DCart
-      USE moduleBoundary
       IMPLICIT NONE
 
       CLASS(meshGeneric), INTENT(inout):: self

src/modules/mesh/makefile

@@ -1,4 +1,4 @@
-all: moduleMesh.o moduleMeshBoundary.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o
+all: moduleMesh.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o
 
 3DCart.o: moduleMesh.o
 	$(MAKE) -C 3DCart all
@@ -20,9 +20,8 @@ all: moduleMesh.o moduleMeshBoundary.o inout.o 3DCart.o 2DCyl.o 2DCart.o 1DRad.o
 
 moduleMesh.o: moduleMesh.f90
 	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
-
+	$(FC) $(FCFLAGS) -c moduleMesh@elements.f90 -o $(OBJDIR)/$@
-moduleMeshBoundary.o: moduleMesh.o moduleMeshBoundary.f90
+	$(FC) $(FCFLAGS) -c moduleMesh@boundary.f90 -o $(OBJDIR)/$@
-	$(FC) $(FCFLAGS) -c $(subst .o,.f90,$@) -o $(OBJDIR)/$@
 
 inout.o: 3DCart.o 2DCyl.o 2DCart.o 1DRad.o 1DCart.o 0D.o
 	$(MAKE) -C inout all

src/modules/mesh/moduleMesh@boundary.f90

@@ -1,9 +1,102 @@
 !moduleMeshBoundary: Boundary functions for the mesh edges
-MODULE moduleMeshBoundary
+submodule(moduleMesh) boundary
-  USE moduleMesh
-
   CONTAINS
-    SUBROUTINE reflection(edge, part)
+    FUNCTION getBoundaryId(physicalSurface) RESULT(id)
+      IMPLICIT NONE
+
+      INTEGER:: physicalSurface
+      INTEGER:: id
+      INTEGER:: i
+
+      id = 0
+      DO i = 1, nBoundary
+        IF (physicalSurface == boundaries(i)%physicalSurface) id = boundaries(i)%n
+
+      END DO
+
+    END FUNCTION getBoundaryId
+
+    SUBROUTINE initWallTemperature(boundary, T, c)
+      USE moduleRefParam
+      IMPLICIT NONE
+
+      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
+      REAL(8), INTENT(in):: T, c !Wall temperature and specific heat
+      REAL(8):: vTh
+
+      vTh = DSQRT(c * T) / v_ref
+      boundary = boundaryWallTemperature(vTh = vTh)
+
+    END SUBROUTINE initWallTemperature
+
+    SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
+      USE moduleRefParam
+      USE moduleSpecies
+      USE moduleCaseParam
+      USE moduleConstParam
+      USE moduleErrors
+      IMPLICIT NONE
+
+      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
+      REAL(8), INTENT(in):: me !Electron mass
+      REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties
+      INTEGER, INTENT(in):: ion
+      INTEGER, OPTIONAL, INTENT(in):: electronSecondary
+      REAL(8):: effTime
+      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection
+      REAL(8), INTENT(in):: eThreshold
+
+      ALLOCATE(boundaryIonization:: boundary)
+
+      SELECT TYPE(boundary)
+      TYPE IS(boundaryIonization)
+        boundary%m0 = m0 / m_ref
+        boundary%n0 = n0 * Vol_ref
+        boundary%v0 = v0 / v_ref
+        boundary%vTh = DSQRT(kb*T0/m0)/v_ref
+        boundary%species => species(ion)%obj
+        IF (PRESENT(electronSecondary)) THEN
+          SELECT TYPE(sp => species(electronSecondary)%obj)
+          TYPE IS(speciesCharged)
+            boundary%electronSecondary => sp
+
+          CLASS DEFAULT
+            CALL criticalError("Species " // sp%name // " chosen for " // &
+                               "secondary electron is not a charged species", 'initIonization')
+
+          END SELECT
+
+        ELSE
+          boundary%electronSecondary => NULL()
+
+        END IF
+        boundary%effectiveTime = effTime / ti_ref
+        CALL boundary%crossSection%init(crossSection)
+        CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
+        boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2)
+        boundary%deltaV = DSQRT(boundary%eThreshold/me)
+
+      END SELECT
+
+    END SUBROUTINE initIonization
+
+    subroutine initQuasiNeutrality(boundary)
+      implicit none
+
+      class(boundaryGeneric), allocatable, intent(out):: boundary
+      integer:: e, et
+
+      allocate(boundaryQuasiNeutrality:: boundary)
+
+      select type(boundary)
+      type is(boundaryQuasiNeutrality)
+        boundary%alpha = 0.d0
+
+      end select
+
+    end subroutine initQuasiNeutrality
+
+    module SUBROUTINE reflection(edge, part)
       USE moduleCaseParam
       USE moduleSpecies
       IMPLICIT NONE
@@ -224,18 +317,23 @@ MODULE moduleMeshBoundary
       class(meshEdge), intent(inout):: edge
       class(particle), intent(inout):: part
       real(8), allocatable:: density(:)
-      integer:: nNodes
+      class(meshCell), pointer:: cell
-      integer, allocatable:: nodes(:)
+      real(8):: EF_dir
+      real(8):: alpha
 
 
       select type(bound => edge%boundary%bTypes(part%species%n)%obj)
       type is(boundaryQuasiNeutrality)
-        nNodes = edge%nNodes
-        nodes  = edge%getNodes(nNodes)
 
-        bound%alpha = 1.0d-1
+        if (associated(edge%e1)) then
+          cell => edge%e1
 
-        if (random() <= bound%alpha) then
+        else
+          cell => edge%e2
+
+        end if
+        
+        if (random() <= alpha) then
           call reflection(edge, part)
 
         else
@@ -248,7 +346,7 @@ MODULE moduleMeshBoundary
     end subroutine quasiNeutrality
 
     !Points the boundary function to specific type
-    SUBROUTINE pointBoundaryFunction(edge, s)
+    module SUBROUTINE pointBoundaryFunction(edge, s)
       USE moduleErrors
       IMPLICIT NONE
 
@@ -284,4 +382,4 @@ MODULE moduleMeshBoundary
 
     END SUBROUTINE pointBoundaryFunction
 
-END MODULE moduleMeshBoundary
+end submodule boundary

src/modules/mesh/moduleMesh@common.f90

@@ -0,0 +1,3 @@
+module moduleMeshCommon
+
+end module moduleMeshCommon

src/modules/mesh/moduleMesh@elements.f90

@@ -0,0 +1,795 @@
+submodule(moduleMesh) elements
+  CONTAINS
+    !Reset the output of node
+    PURE module SUBROUTINE resetOutput(self)
+      USE moduleSpecies
+      USE moduleOutput
+      IMPLICIT NONE
+
+      CLASS(meshNode), INTENT(inout):: self
+      INTEGER:: k
+
+      DO k = 1, nSpecies
+        self%output(k)%den     = 0.D0
+        self%output(k)%mom     = 0.D0
+        self%output(k)%tensorS = 0.D0
+
+      END DO
+
+    END SUBROUTINE resetOutput
+
+    !Constructs the global K matrix
+    PURE module SUBROUTINE constructGlobalK(self)
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(inout):: self
+      INTEGER:: c
+      INTEGER, ALLOCATABLE:: n(:)
+      REAL(8), ALLOCATABLE:: localK(:,:)
+      INTEGER:: i, j
+
+      DO c = 1, self%numCells
+        associate(nNodes => self%cells(c)%obj%nNodes)
+          ALLOCATE(n(1:nNodes))
+          ALLOCATE(localK(1:nNodes, 1:nNodes))
+          n      = self%cells(c)%obj%getNodes(nNodes)
+          localK = self%cells(c)%obj%elemK(nNodes)
+
+          DO i = 1, nNodes
+            DO j = 1, nNodes
+              self%K(n(i), n(j)) = self%K(n(i), n(j)) + localK(i, j)
+
+            END DO
+
+          END DO
+
+          DEALLOCATE(n, localK)
+
+        end associate
+
+      END DO
+
+    END SUBROUTINE constructGlobalK
+
+    ! Gather the value of valNodes at position Xi of an edge
+    pure module function gatherF_edge_scalar(self, Xi, nNodes, valNodes) RESULT(f)
+      implicit none
+
+      class(meshEdge), intent(in):: self
+      real(8), intent(in):: Xi(1:3)
+      integer, intent(in):: nNodes
+      real(8), intent(in):: valNodes(1:nNodes)
+      real(8):: f
+      real(8):: fPsi(1:nNodes)
+
+      fPsi = self%fPsi(Xi, nNodes)
+      f = dot_product(fPsi, valNodes)
+
+    end function gatherF_edge_scalar
+
+    !Gather the value of valNodes (scalar) at position Xi
+    PURE module FUNCTION gatherF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(f)
+      IMPLICIT NONE
+
+      CLASS(meshCell), INTENT(in):: self
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8), INTENT(in):: valNodes(1:nNodes)
+      REAL(8):: f
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi = self%fPsi(Xi, nNodes)
+      f = DOT_PRODUCT(fPsi, valNodes)
+
+    END FUNCTION gatherF_cell_scalar
+
+    !Gather the value of valNodes (array) at position Xi
+    PURE module FUNCTION gatherF_cell_array(self, Xi, nNodes, valNodes) RESULT(f)
+      IMPLICIT NONE
+
+      CLASS(meshCell), INTENT(in):: self
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8), INTENT(in):: valNodes(1:nNodes, 1:3)
+      REAL(8):: f(1:3)
+      REAL(8):: fPsi(1:nNodes)
+
+      fPsi = self%fPsi(Xi, nNodes)
+      f = MATMUL(fPsi, valNodes)
+
+    END FUNCTION gatherF_cell_array
+
+    !Gather the spatial derivative of valNodes (scalar) at position Xi
+    PURE module FUNCTION gatherDF_cell_scalar(self, Xi, nNodes, valNodes) RESULT(df)
+      IMPLICIT NONE
+
+      CLASS(meshCell), INTENT(in):: self
+      REAL(8), INTENT(in):: Xi(1:3)
+      INTEGER, INTENT(in):: nNodes
+      REAL(8), INTENT(in):: valNodes(1:nNodes)
+      REAL(8):: df(1:3)
+      REAL(8):: dPsi(1:3, 1:nNodes)
+      REAL(8):: pDer(1:3,1:3)
+      REAL(8):: dPsiR(1:3, 1:nNodes)
+      REAL(8):: invJ(1:3, 1:3), detJ
+
+      dPsi  = self%dPsi(Xi, nNodes)
+      pDer  = self%partialDer(nNodes, dPsi)
+      detJ  = self%detJac(pDer)
+      invJ  = self%invJac(pDer)
+      dPsiR =  MATMUL(invJ, dPsi)/detJ
+      df    = (/ DOT_PRODUCT(dPsiR(1,:), valNodes), &
+                 DOT_PRODUCT(dPsiR(2,:), valNodes), &
+                 DOT_PRODUCT(dPsiR(3,:), valNodes) /)
+
+    END FUNCTION gatherDF_cell_scalar
+
+    !Scatters particle properties into cell nodes
+    module SUBROUTINE scatter(self, nNodes, part)
+      USE moduleMath
+      USE moduleSpecies
+      USE OMP_LIB
+      IMPLICIT NONE
+
+      CLASS(meshCell), INTENT(inout):: self
+      INTEGER, INTENT(in):: nNodes
+      CLASS(particle), INTENT(in):: part
+      REAL(8):: fPsi(1:nNodes)
+      INTEGER:: cellNodes(1:nNodes)
+      REAL(8):: tensorS(1:3, 1:3)
+      INTEGER:: sp
+      INTEGER:: i
+      CLASS(meshNode), POINTER:: node
+      REAL(8):: pFraction !Particle fraction
+
+      cellNodes = self%getNodes(nNodes)
+      fPsi = self%fPsi(part%Xi, nNodes)
+
+      tensorS = outerProduct(part%v, part%v)
+
+      sp = part%species%n
+
+      DO i = 1, nNodes
+        node => mesh%nodes(cellNodes(i))%obj
+        pFraction = fPsi(i)*part%weight
+        CALL OMP_SET_LOCK(node%lock)
+        node%output(sp)%den          = node%output(sp)%den          + pFraction
+        node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + pFraction*part%v(:)
+        node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS
+        CALL OMP_UNSET_LOCK(node%lock)
+
+      END DO
+
+    END SUBROUTINE scatter
+
+    !Find next cell for particle
+    RECURSIVE module SUBROUTINE findCell(self, part, oldCell)
+      USE moduleSpecies
+      USE moduleErrors
+      USE OMP_LIB
+      IMPLICIT NONE
+
+      CLASS(meshCell), INTENT(inout):: self
+      CLASS(particle), INTENT(inout), TARGET:: part
+      CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell
+      REAL(8):: Xi(1:3)
+      CLASS(meshElement), POINTER:: neighbourElement
+      INTEGER:: sp
+
+      Xi = self%phy2log(part%r)
+      !Checks if particle is inside 'self' cell
+      IF (self%inside(Xi)) THEN
+        part%cell = self%n
+        part%Xi   = Xi
+        part%n_in = .TRUE.
+        !Assign particle to listPart_in
+        IF (listInCells) THEN
+          CALL OMP_SET_LOCK(self%lock)
+          sp = part%species%n
+          CALL self%listPart_in(sp)%add(part)
+          self%totalWeight(sp) = self%totalWeight(sp) + part%weight
+          CALL OMP_UNSET_LOCK(self%lock)
+
+        END IF
+
+      ELSE
+        !If not, searches for a neighbour and repeats the process.
+        CALL self%neighbourElement(Xi, neighbourElement)
+        !Defines the next step
+        SELECT TYPE(neighbourElement)
+        CLASS IS(meshCell)
+          !Particle moved to new cell, repeat find procedure
+          CALL neighbourElement%findCell(part, self)
+
+        CLASS IS (meshEdge)
+          !Particle encountered a surface, apply boundary
+          CALL neighbourElement%fBoundary(part%species%n)%apply(neighbourElement,part)
+
+          !If particle is still inside the domain, call findCell
+          IF (part%n_in) THEN
+            IF(PRESENT(oldCell)) THEN
+              CALL self%findCell(part, oldCell)
+
+            ELSE
+              CALL self%findCell(part)
+
+            END IF
+          END IF
+
+        CLASS DEFAULT
+          WRITE (*, "(A, I6)") "Element = ", self%n
+          CALL criticalError("No connectivity found for element", "findCell")
+
+        END SELECT
+
+      END IF
+
+    END SUBROUTINE findCell
+
+    !If Coll and Particle are the same, simply copy the part%cell into part%cellColl
+    module SUBROUTINE findCellSameMesh(part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+
+      part%cellColl = part%cell
+
+    END SUBROUTINE findCellSameMesh
+
+    !TODO: try to combine this with the findCell for a regular mesh
+    !Find the volume in which particle reside in the mesh for collisions
+    !No boundary interaction taken into account
+    module SUBROUTINE findCellCollMesh(part)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      TYPE(particle), INTENT(inout):: part
+      LOGICAL:: found
+      CLASS(meshCell), POINTER:: cell
+      REAL(8), DIMENSION(1:3):: Xi
+      CLASS(meshElement), POINTER:: neighbourElement
+      INTEGER:: sp
+
+      found = .FALSE.
+
+      cell => meshColl%cells(part%cellColl)%obj
+      DO WHILE(.NOT. found)
+        Xi = cell%phy2log(part%r)
+        IF (cell%inside(Xi)) THEN
+          part%cellColl = cell%n
+          IF (listInCells) THEN
+            CALL OMP_SET_LOCK(cell%lock)
+            sp = part%species%n
+            CALL cell%listPart_in(sp)%add(part)
+            cell%totalWeight(sp) = cell%totalWeight(sp) + part%weight
+            CALL OMP_UNSET_LOCK(cell%lock)
+
+          END IF
+          found = .TRUE.
+
+        ELSE
+          CALL cell%neighbourElement(Xi, neighbourElement)
+          SELECT TYPE(neighbourElement)
+          CLASS IS(meshCell)
+            !Try next element
+            cell => neighbourElement
+
+          CLASS DEFAULT
+            !Should never happend, but just in case, stops loops
+            found = .TRUE.
+
+          END SELECT
+
+        END IF
+
+      END DO
+
+    END SUBROUTINE findCellCollMesh
+
+    !Returns index of volume associated to a position (if any)
+    !If no voulme is found, returns 0
+    !WARNING: This function is slow and should only be used in initialization phase
+    module FUNCTION findCellBrute(self, r) RESULT(nVol)
+      USE moduleSpecies
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(in):: self
+      REAL(8), DIMENSION(1:3), INTENT(in):: r
+      INTEGER:: nVol
+      INTEGER:: e
+      REAL(8), DIMENSION(1:3):: Xi
+
+      !Inits RESULT
+      nVol = 0
+
+      DO e = 1, self%numCells
+        Xi = self%cells(e)%obj%phy2log(r)
+        IF(self%cells(e)%obj%inside(Xi)) THEN
+          nVol = self%cells(e)%obj%n
+          EXIT
+
+        END IF
+
+      END DO
+
+    END FUNCTION findCellBrute
+
+    !Computes collisions in element
+    module SUBROUTINE doCollisions(self)
+      USE moduleCollisions
+      USE moduleSpecies
+      USE moduleList
+      use moduleRefParam
+      USE moduleRandom
+      USE moduleOutput
+      USE moduleMath
+      USE moduleCaseParam, ONLY: timeStep
+      IMPLICIT NONE
+
+      CLASS(meshGeneric), INTENT(inout), TARGET:: self
+      INTEGER:: e
+      CLASS(meshCell), POINTER:: cell
+      INTEGER:: k, i, j
+      INTEGER:: nPart_i, nPart_j, nPart!Number of particles inside the cell
+      REAL(8):: pMax !Maximum probability of collision
+      INTEGER:: nColl
+      TYPE(pointerArray), ALLOCATABLE:: partTemp_i(:), partTemp_j(:)
+      TYPE(particle), POINTER:: part_i, part_j
+      INTEGER:: n, c
+      REAL(8):: vRel, rMass, eRel
+      REAL(8):: sigmaVrelTotal
+      REAL(8), ALLOCATABLE:: sigmaVrel(:), probabilityColl(:)
+      REAL(8):: rnd_real !Random number for collision
+      INTEGER:: rnd_int  !Random number for collision
+
+      IF (MOD(timeStep, everyColl) == 0) THEN
+        !Collisions need to be performed in this iteration
+        !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(part_i, part_j, partTemp_i, partTemp_j)
+        DO e=1, self%numCells
+
+          cell => self%cells(e)%obj
+
+          !TODO: Simplify this, to many sublevels
+          !Iterate over the number of pairs
+          DO k = 1, nCollPairs
+            !Reset tally of collisions
+            IF (collOutput) THEN
+              cell%tallyColl(k)%tally = 0
+
+            END IF
+
+            IF (interactionMatrix(k)%amount > 0) THEN
+              !Select the species for the collision pair
+              i = interactionMatrix(k)%sp_i%n
+              j = interactionMatrix(k)%sp_j%n
+
+              !Number of particles per species in the collision pair
+              nPart_i = cell%listPart_in(i)%amount
+              nPart_j = cell%listPart_in(j)%amount
+
+              IF (nPart_i > 0 .AND. nPart_j > 0) THEN
+                !Total number of particles for the collision pair
+                nPart = nPart_i + nPart_j
+
+                !Resets the number of collisions in the cell
+                nColl = 0
+
+                !Probability of collision for pair i-j
+                pMax = (cell%totalWeight(i) + cell%totalWeight(j))*cell%sigmaVrelMax(k)*tauColl/cell%volume
+
+                !Number of collisions in the cell
+                nColl = NINT(REAL(nPart)*pMax*0.5D0)
+
+                !Converts the list of particles to an array for easy access
+                IF (nColl > 0) THEN
+                  partTemp_i = cell%listPart_in(i)%convert2Array()
+                  partTemp_j = cell%listPart_in(j)%convert2Array()
+
+                END IF
+
+                DO n = 1, nColl
+                  !Select random particles
+                  part_i => NULL()
+                  part_j => NULL()
+                  rnd_int = random(1, nPart_i)
+                  part_i => partTemp_i(rnd_int)%part
+                  rnd_int = random(1, nPart_j)
+                  part_j => partTemp_j(rnd_int)%part
+                  !If they are the same particle, skip
+                  !TODO: Maybe try to improve this
+                  IF (ASSOCIATED(part_i, part_j)) THEN
+                    CYCLE
+
+                  END IF
+
+                  !If particles do not belong to the species, skip collision
+                  !This can happen, for example, if particle has been previously ionized or removed
+                  !TODO: Try to find a way to not lose these collisions. Maybe check new 'k' and use that for the collision?
+                  IF (part_i%species%n /= i .OR. &
+                      part_j%species%n /= j) THEN
+                      CYCLE
+
+                  END IF
+                  !Obtain the cross sections for the different processes
+                  !TODO: From here it might be a procedure in interactionMatrix
+                  vRel = NORM2(part_i%v-part_j%v)
+                  rMass = reducedMass(part_i%weight*part_i%species%m, part_j%weight*part_j%species%m)
+                  eRel = rMass*vRel**2
+                  CALL interactionMatrix(k)%getSigmaVrel(vRel, eRel, sigmaVrelTotal, sigmaVrel)
+
+                  !Update maximum sigma*v_rel
+                  IF (sigmaVrelTotal > cell%sigmaVrelMax(k)) THEN
+                    cell%sigmaVrelMax(k) = sigmaVrelTotal
+
+                  END IF
+
+                  ALLOCATE(probabilityColl(0:interactionMatrix(k)%amount))
+                  probabilityColl = 0.0
+                  DO c = 1, interactionMatrix(k)%amount
+                    probabilityColl(c) = sigmaVrel(c)/cell%sigmaVrelMax(k) + SUM(probabilityColl(0:c-1))
+
+                  END DO
+
+                  !Selects random number between 0 and 1
+                  rnd_real = random()
+
+                  !If the random number is below the total probability of collision, collide particles
+                  IF (rnd_real < sigmaVrelTotal / cell%sigmaVrelMax(k)) THEN
+
+                    !Loop over collisions
+                    DO c = 1, interactionMatrix(k)%amount
+                      IF (rnd_real <= probabilityColl(c)) THEN
+                        !$OMP CRITICAL
+                        CALL interactionMatrix(k)%collisions(c)%obj%collide(part_i, part_j, vRel)
+                        !$OMP END CRITICAL
+
+                        !If collisions are gonna be output, count the collision
+                        IF (collOutput) THEN
+                          cell%tallyColl(k)%tally(c) = cell%tallyColl(k)%tally(c) + 1
+
+                        END IF
+
+                        !A collision has ocurred, exit the loop
+                        EXIT
+
+                      END IF
+
+                    END DO
+
+                  END IF
+
+                  !Deallocate arrays for next collision
+                  DEALLOCATE(sigmaVrel, probabilityColl)
+
+                !End loop collisions in cell
+                END DO
+
+              END IF
+
+            END IF
+
+          !End loop collision pairs
+          END DO
+
+        !End loop volumes
+        END DO
+        !$OMP END DO
+
+      END IF
+
+    END SUBROUTINE doCollisions
+
+    module SUBROUTINE doCoulomb(self)
+      USE moduleCoulomb
+      USE moduleRandom
+      USE moduleOutput
+      USE moduleList
+      USE moduleMath
+      USE moduleRefParam
+      USE moduleConstParam
+      IMPLICIT NONE
+
+      CLASS(meshParticles), INTENT(in), TARGET:: self
+      CLASS(meshCell), POINTER:: cell
+      TYPE(interactionsCoulomb):: pair
+      INTEGER:: e
+      INTEGER:: k
+      INTEGER:: i, j
+      INTEGER:: n
+      INTEGER:: p
+      TYPE(lNode), POINTER:: partTemp
+      INTEGER(8), ALLOCATABLE:: cellNodes(:)
+      CLASS(meshNode), POINTER:: node
+      TYPE(outputFormat):: output
+      REAL(8), ALLOCATABLE:: densityNodes(:), velocityNodes(:,:), temperatureNodes(:) !values in node
+      REAL(8):: density, velocity(1:3), temperature!values at particle position
+      REAL(8):: C(1:3), C_per, W(1:3) !relative velocity and velocity in the relative frame of reference
+      REAL(8):: l, lW, l2
+      REAL(8):: GlW, HlW
+      REAL(8):: normC
+      REAL(8):: cosThe, sinThe
+      REAL(8):: cosPhi, sinPhi
+      REAL(8):: rotation(1:3,1:3) !Rotation matrix to go back to laboratory frame
+      REAL(8):: A, AW
+      REAL(8):: deltaW_par, deltaW_par_square, deltaW_per_square !Increments of W
+      REAL(8):: theta_per !Random angle for perpendicular direction
+      REAL(8):: eps = 1.D-12
+      REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ij, p_ij
+      REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ij
+      REAL(8):: massSum_ij
+      REAL(8), ALLOCATABLE, DIMENSION(:,:):: deltaV_ji, p_ji
+      REAL(8), ALLOCATABLE, DIMENSION(:):: mass_ji
+      REAL(8):: massSum_ji
+      REAL(8):: alpha_num, alpha_den, alpha, beta(1:3)
+
+
+      !$OMP DO SCHEDULE(DYNAMIC) PRIVATE(partTemp)
+      DO e = 1, self%numCells
+        cell => self%cells(e)%obj
+        cellNodes = cell%getNodes(cell%nNodes)
+
+        ALLOCATE(densityNodes(1:cell%nNodes),       &
+                 velocityNodes(1:cell%nNodes, 1:3), &
+                 temperatureNodes(1:cell%nNodes))
+
+        DO k=1, nCoulombPairs
+          pair = coulombMatrix(k)
+          i = pair%sp_i%n
+          j = pair%sp_j%n
+
+          !Do scattering of particles from species_i due to species j
+          !Compute background properties of species_j
+          DO n = 1, cell%nNodes
+            node => self%nodes(cellNodes(n))%obj
+            CALL calculateOutput(node%output(j), output, node%v, pair%sp_j)
+            densityNodes(n)      = output%density/n_ref
+            velocityNodes(n,1:3) = output%velocity(1:3)/v_ref
+            temperatureNodes(n)  = output%temperature/T_ref
+
+          END DO
+
+          ALLOCATE(deltaV_ij(1:cell%listPart_in(i)%amount, 1:3))
+          ALLOCATE(p_ij(1:cell%listPart_in(i)%amount, 1:3))
+          ALLOCATE(mass_ij(1:cell%listPart_in(i)%amount))
+          deltaV_ij = 0.D0
+          p_ij      = 0.D0
+          mass_ij   = 0.D0
+          !Loop over particles of species_i
+          partTemp => cell%listPart_in(i)%head
+          p = 1
+          DO WHILE(ASSOCIATED(partTemp))
+            density     = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes)
+            velocity    = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes)
+            temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes)
+
+            !If cell temperature is too low, skip particle to avoid division by zero
+            IF (temperature>eps) THEN
+              l2 = pair%l2_j/temperature
+              l  = SQRT(l2)
+
+            ELSE
+              partTemp => partTemp%next
+
+              CYCLE
+
+            END IF
+
+            A  = pair%A_i*density
+
+            C = partTemp%part%v - velocity
+            normC = NORM2(C)
+
+            !C_3 = z; C_1, C2 = x, y (per)
+            C_per = NORM2(C(1:2))
+            cosPhi = C(1)  / C_per
+            sinPhi = C(2)  / C_per
+            cosThe = C(3)  / normC
+            sinThe = C_per / normC
+
+            !Rotation matrix to go from W to C
+            rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe,    & !First column
+                                        -sinPhi,        cosPhi,    0.D0,    & !Second column
+                                  sinThe*cosPhi, sinThe*sinPhi,  cosThe /), & !Third column
+                                  (/ 3, 3 /))
+
+            !W at start is = (0, 0, normC), so normW = normC
+            lW = l * normC
+            GlW = G(lW)
+            HlW = H(lW)
+            AW = A / normC
+
+            !Calculate changes in W due to collision process
+            deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin
+            deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian()
+            deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian()
+
+            !Random angle to distribute perpendicular change in velocity
+            theta_per = PI2*random()
+
+            !Change W
+            W(1) =         deltaW_per_square * COS(theta_per) 
+            W(2) =         deltaW_per_square * SIN(theta_per)
+            W(3) = normC + deltaW_par + deltaW_par_square
+
+            !Compute changes in velocity for each particle
+            deltaV_ij(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
+            mass_ij(p)       = pair%sp_i%m*partTemp%part%weight
+            p_ij(p,1:3)      = mass_ij(p)*partTemp%part%v
+
+            !Move to the next particle in the list
+            partTemp => partTemp%next
+            p = p + 1
+
+          END DO
+
+          !Do corresponding collisions
+          IF (i /= j) THEN
+            !Do scattering of particles from species_j due to species i
+            !Compute background properties of species_i
+            DO n = 1, cell%nNodes
+              node => self%nodes(cellNodes(n))%obj
+              CALL calculateOutput(node%output(i), output, node%v, pair%sp_i)
+              densityNodes(n)      = output%density/n_ref
+              velocityNodes(n,1:3) = output%velocity(1:3)/v_ref
+              temperatureNodes(n)  = output%temperature/T_ref
+
+            END DO
+
+            ALLOCATE(deltaV_ji(1:cell%listPart_in(j)%amount, 1:3))
+            ALLOCATE(p_ji(1:cell%listPart_in(j)%amount, 1:3))
+            ALLOCATE(mass_ji(1:cell%listPart_in(j)%amount))
+            deltaV_ji = 0.D0
+            p_ji      = 0.D0
+            mass_ji   = 0.D0
+            !Loop over particles of species_j
+            partTemp => cell%listPart_in(j)%head
+            p = 1
+            DO WHILE(ASSOCIATED(partTemp))
+              density     = cell%gatherF(partTemp%part%Xi, cell%nNodes, densityNodes)
+              velocity    = cell%gatherF(partTemp%part%Xi, cell%nNodes, velocityNodes)
+              temperature = cell%gatherF(partTemp%part%Xi, cell%nNodes, temperatureNodes)
+
+              !If cell temperature is too low, skip particle to avoid division by zero
+              IF (temperature>eps) THEN
+                l2 = pair%l2_i/temperature
+                l  = SQRT(l2)
+
+              ELSE
+                partTemp => partTemp%next
+
+                CYCLE
+
+              END IF
+              A  = pair%A_j*density
+
+              C = partTemp%part%v - velocity
+              normC = NORM2(C)
+
+              !C_3 = z; C_1, C2 = x, y (per)
+              C_per = NORM2(C(1:2))
+              cosPhi = C(1)  / C_per
+              sinPhi = C(2)  / C_per
+              cosThe = C(3)  / normC
+              sinThe = C_per / normC
+
+              !Rotation matrix to go from W to C
+              rotation = RESHAPE((/ cosThe*cosPhi, cosThe*sinPhi, -sinThe,    & !First column
+                                          -sinPhi,        cosPhi,    0.D0,    & !Second column
+                                    sinThe*cosPhi, sinThe*sinPhi,  cosThe /), & !Third column
+                                    (/ 3, 3 /))
+
+              !W at start is = (0, 0, normC), so normW = normC
+              lW = l * normC
+              GlW = G(lW)
+              HlW = H(lW)
+              AW = A / normC
+
+              !Calculate changes in W due to collision process
+              deltaW_par = - A * pair%one_plus_massRatio_ij * l2 * GlW * tauMin
+              deltaW_par_square = SQRT(AW * GlW * tauMin)*randomMaxwellian()
+              deltaW_per_square = SQRT(AW * HlW * tauMin)*randomMaxwellian()
+
+              !Random angle to distribute perpendicular change in velocity
+              theta_per = PI2*random()
+
+              !Change W
+              W(1) =         deltaW_per_square * COS(theta_per) 
+              W(2) =         deltaW_per_square * SIN(theta_per)
+              W(3) = normC + deltaW_par + deltaW_par_square
+
+              !Compute changes in velocity for each particle
+              deltaV_ji(p,1:3) = MATMUL(rotation, W) + velocity - partTemp%part%v
+              mass_ji(p)       = pair%sp_j%m*partTemp%part%weight
+              p_ji(p,1:3)      = mass_ji(p)*partTemp%part%v
+
+              !Move to the next particle in the list
+              partTemp => partTemp%next
+              p = p + 1
+
+            END DO
+
+          END IF
+
+          !Calculate correction
+          !Total mass
+          massSum_ij = SUM(mass_ij)
+          massSum_ji = 0.D0
+
+          !Beta
+          beta = 0.D0
+          DO p = 1, cell%listPart_in(i)%amount
+            beta = beta + mass_ij(p) * deltaV_ij(p,1:3)
+
+          END DO
+          
+          IF (i /= j) THEN
+            massSum_ji = SUM(mass_ji)
+            DO p = 1, cell%listPart_in(j)%amount
+              beta = beta + mass_ji(p) * deltaV_ji(p,1:3)
+
+            END DO
+
+          END IF
+
+          beta = beta / (massSum_ij + massSum_ji)
+
+          !Alpha
+          alpha_num = 0.D0
+          alpha_den = 0.D0
+          DO p =1, cell%listPart_in(i)%amount
+            alpha_num = alpha_num + DOT_PRODUCT(p_ij(p,1:3), deltav_ij(p,1:3) - beta(1:3))
+            alpha_den = alpha_den + mass_ij(p) * NORM2(deltav_ij(p,1:3) - beta(1:3))**2
+
+          END DO
+
+          IF (i /= j) THEN
+            DO p = 1, cell%listPart_in(j)%amount
+              alpha_num = alpha_num + DOT_PRODUCT(p_ji(p,1:3), deltav_ji(p,1:3) - beta(1:3))
+              alpha_den = alpha_den + mass_ji(p) * NORM2(deltav_ji(p,1:3) - beta(1:3))**2
+
+            END DO
+
+          END IF
+
+          alpha = -2.D0*alpha_num / alpha_den
+
+          !Apply correction to particles velocity
+          partTemp => cell%listPart_in(i)%head
+          p = 1
+          DO WHILE(ASSOCIATED(partTemp))
+            partTemp%part%v = partTemp%part%v + alpha * (deltaV_ij(p,1:3) - beta(1:3))
+            partTemp => partTemp%next
+            p = p + 1
+
+          END DO
+
+          IF (i /= j) THEN
+            partTemp => cell%listPart_in(j)%head
+            p = 1
+            DO WHILE(ASSOCIATED(partTemp))
+              partTemp%part%v = partTemp%part%v + alpha * (deltaV_ji(p,1:3) - beta(1:3))
+              partTemp => partTemp%next
+              p = p + 1
+
+            END DO
+
+          END IF
+
+          DEALLOCATE(deltaV_ij, p_ij, mass_ij)
+
+          IF (i /= j) THEN
+            DEALLOCATE(deltaV_ji, p_ji, mass_ji)
+
+          END IF
+
+        END DO
+
+        DEALLOCATE(densityNodes, velocityNodes, temperatureNodes, cellNodes)
+
+      END DO
+      !$OMP END DO
+
+    END SUBROUTINE doCoulomb
+
+end submodule elements

src/modules/moduleBoundary.f90

@@ -1,178 +0,0 @@
-MODULE moduleBoundary
-  USE moduleTable
-  USE moduleSpecies
-
-  !Generic type for boundaries
-  TYPE, PUBLIC:: boundaryGeneric
-    CONTAINS
-
-  END TYPE boundaryGeneric
-
-  !Reflecting boundary
-  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryReflection
-    CONTAINS
-
-  END TYPE boundaryReflection
-
-  !Absorption boundary
-  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAbsorption
-    CONTAINS
-
-  END TYPE boundaryAbsorption
-
-  !Transparent boundary
-  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryTransparent
-    CONTAINS
-
-  END TYPE boundaryTransparent
-
-  !Symmetry axis
-  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryAxis
-    CONTAINS
-
-  END TYPE boundaryAxis
-
-  !Wall Temperature boundary
-  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryWallTemperature
-    !Thermal velocity of the wall: square root(Wall temperature X specific heat)
-    REAL(8):: vTh
-    CONTAINS
-
-  END TYPE boundaryWallTemperature
-
-  !Ionization boundary
-  TYPE, PUBLIC, EXTENDS(boundaryGeneric):: boundaryIonization
-    REAL(8):: m0, n0, v0(1:3), vTh !Properties of background neutrals.
-    CLASS(speciesGeneric), POINTER:: species !Ion species
-    CLASS(speciesCharged), POINTER:: electronSecondary !Pointer to species considerer as secondary electron
-    TYPE(table1D):: crossSection
-    REAL(8):: effectiveTime
-    REAL(8):: eThreshold
-    REAL(8):: deltaV
-    CONTAINS
-
-  END TYPE boundaryIonization
-
-  ! Ensures quasi-neutrality by changing the reflection coefficient
-  type, public, extends(boundaryGeneric):: boundaryQuasiNeutrality
-    real(8):: alpha ! Reflection parameter
-
-  end type boundaryQuasiNeutrality
-  !Wrapper for boundary types (one per species)
-  TYPE:: bTypesCont
-    CLASS(boundaryGeneric), ALLOCATABLE:: obj
-    CONTAINS
-
-  END TYPE bTypesCont
-
-  !Wrapper for boundary conditions
-  TYPE:: boundaryCont
-    INTEGER:: n = 0
-    CHARACTER(:), ALLOCATABLE:: name 
-    INTEGER:: physicalSurface = 0 !Physical surface as defined in the mesh file
-    CLASS(bTypesCont), ALLOCATABLE:: bTypes(:) !Array for boundary per species
-    CONTAINS
-
-  END TYPE boundaryCont
-
-  !Number of boundaries
-  INTEGER:: nBoundary = 0
-  !Array for boundaries
-  TYPE(boundaryCont), ALLOCATABLE, TARGET:: boundary(:)
-
-  CONTAINS
-    FUNCTION getBoundaryId(physicalSurface) RESULT(id)
-      IMPLICIT NONE
-
-      INTEGER:: physicalSurface
-      INTEGER:: id
-      INTEGER:: i
-
-      id = 0
-      DO i = 1, nBoundary
-        IF (physicalSurface == boundary(i)%physicalSurface) id = boundary(i)%n
-
-      END DO
-
-    END FUNCTION getBoundaryId
-
-    SUBROUTINE initWallTemperature(boundary, T, c)
-      USE moduleRefParam
-      IMPLICIT NONE
-
-      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
-      REAL(8), INTENT(in):: T, c !Wall temperature and specific heat
-      REAL(8):: vTh
-
-      vTh = DSQRT(c * T) / v_ref
-      boundary = boundaryWallTemperature(vTh = vTh)
-
-    END SUBROUTINE initWallTemperature
-
-    SUBROUTINE initIonization(boundary, me, m0, n0, v0, T0, ion, effTime, crossSection, eThreshold, electronSecondary)
-      USE moduleRefParam
-      USE moduleSpecies
-      USE moduleCaseParam
-      USE moduleConstParam
-      USE moduleErrors
-      IMPLICIT NONE
-
-      CLASS(boundaryGeneric), ALLOCATABLE, INTENT(out):: boundary
-      REAL(8), INTENT(in):: me !Electron mass
-      REAL(8), INTENT(in):: m0, n0, v0(1:3), T0 !Neutral properties
-      INTEGER, INTENT(in):: ion
-      INTEGER, OPTIONAL, INTENT(in):: electronSecondary
-      REAL(8):: effTime
-      CHARACTER(:), ALLOCATABLE, INTENT(in):: crossSection
-      REAL(8), INTENT(in):: eThreshold
-
-      ALLOCATE(boundaryIonization:: boundary)
-
-      SELECT TYPE(boundary)
-      TYPE IS(boundaryIonization)
-        boundary%m0 = m0 / m_ref
-        boundary%n0 = n0 * Vol_ref
-        boundary%v0 = v0 / v_ref
-        boundary%vTh = DSQRT(kb*T0/m0)/v_ref
-        boundary%species => species(ion)%obj
-        IF (PRESENT(electronSecondary)) THEN
-          SELECT TYPE(sp => species(electronSecondary)%obj)
-          TYPE IS(speciesCharged)
-            boundary%electronSecondary => sp
-
-          CLASS DEFAULT
-            CALL criticalError("Species " // sp%name // " chosen for " // &
-                               "secondary electron is not a charged species", 'initIonization')
-
-          END SELECT
-
-        ELSE
-          boundary%electronSecondary => NULL()
-
-        END IF
-        boundary%effectiveTime = effTime / ti_ref
-        CALL boundary%crossSection%init(crossSection)
-        CALL boundary%crossSection%convert(eV2J/(m_ref*v_ref**2), 1.D0/L_ref**2)
-        boundary%eThreshold = eThreshold*eV2J/(m_ref*v_ref**2)
-        boundary%deltaV = DSQRT(boundary%eThreshold/me)
-
-      END SELECT
-
-    END SUBROUTINE initIonization
-
-    subroutine initQuasiNeutrality(boundary)
-      implicit none
-
-      class(boundaryGeneric), allocatable, intent(out):: boundary
-
-      allocate(boundaryQuasiNeutrality:: boundary)
-
-      select type(boundary)
-      type is(boundaryQuasiNeutrality)
-        boundary%alpha = 0.d0
-
-      end select
-
-    end subroutine initQuasiNeutrality
-
-END MODULE moduleBoundary

src/modules/moduleInject.f90

@@ -324,7 +324,6 @@ MODULE moduleInject
       USE moduleMesh
       USE moduleRandom
       USE moduleErrors
-      use moduleMeshBoundary, only: reflection
       IMPLICIT NONE
 
       CLASS(injectGeneric), INTENT(in):: self