Injection of particles from surfaces direction

Now, if no normal is provided to an injection in the input file, the
velocity direction of the particles is chosen to be the surface normal.

This allows to inject particles from curves, corners... without having
to provide a direction or declaring multiple injections.

doc/user-manual/fpakc_UserManual.tex

@@ -608,7 +608,8 @@ make
           \item \textbf{n}: Real.
                             Array dimension $3$.
                             Direction of injection.
-                            This vector is normalized to $1$.
+                            If no value is provided, the normal of the edge is used as the direction of injection.
+                            This vector is normalized to $1$ by the code.
           \item \textbf{velDist}: Character.
                                   Array dimension $3$.
                                   Type of distribution function used to obtain injected particle velocity:

src/modules/init/moduleInput.f90

@@ -1140,6 +1140,7 @@ MODULE moduleInput
       CHARACTER(:), ALLOCATABLE:: name
       REAL(8):: v
       REAL(8), ALLOCATABLE:: T(:), normal(:)
+      LOGICAL:: fixDirection
       REAL(8):: flow
       CHARACTER(:), ALLOCATABLE:: units
       INTEGER:: physicalSurface
@@ -1160,11 +1161,18 @@ MODULE moduleInput
         CALL config%get(object // '.v',               v,               found)
         CALL config%get(object // '.T',               T,               found)
         CALL config%get(object // '.n',               normal,          found)
+        IF (found) THEN
+          fixDirection = .TRUE.
+        ELSE
+          fixDirection = .FALSE.
+          ALLOCATE(normal(1:3))
+          normal = 0.D0
+        END IF
         CALL config%get(object // '.flow',            flow,            found)
         CALL config%get(object // '.units',           units,           found)
         CALL config%get(object // '.physicalSurface', physicalSurface, found)
 
-        CALL inject(i)%init(i, v, normal, T, flow, units, sp, physicalSurface)
+        CALL inject(i)%init(i, v, normal, fixDirection, T, flow, units, sp, physicalSurface)
 
         CALL readVelDistr(config, inject(i), object)
         
@@ -1240,13 +1248,12 @@ MODULE moduleInput
 
         SELECT CASE(fvType)
         CASE ("Maxwellian")
-          v = inj%vMod*inj%n(i)
           T = inj%T(i)
-          CALL initVelDistMaxwellian(inj%v(i)%obj, v, t, m)
+          CALL initVelDistMaxwellian(inj%v(i)%obj, t, m)
 
         CASE ("Delta")
           v = inj%vMod*inj%n(i)
-          CALL initVelDistDelta(inj%v(i)%obj, v)
+          CALL initVelDistDelta(inj%v(i)%obj)
 
         CASE DEFAULT
           CALL criticalError("No velocity distribution type " // fvType // " defined", 'readVelDistr')

src/modules/moduleInject.f90

@@ -29,7 +29,6 @@ MODULE moduleInject
 
   !Maxwellian distribution function
   TYPE, EXTENDS(velDistGeneric):: velDistMaxwellian
-    REAL(8):: v !Velocity
     REAL(8):: vTh !Thermal Velocity
     CONTAINS
       PROCEDURE, PASS:: randomVel => randomVelMaxwellian
@@ -38,7 +37,6 @@ MODULE moduleInject
 
   !Dirac's delta distribution function
   TYPE, EXTENDS(velDistGeneric):: velDistDelta
-    REAL(8):: v !Velocity
     CONTAINS
       PROCEDURE, PASS:: randomVel => randomVelDelta
 
@@ -51,6 +49,7 @@ MODULE moduleInject
     REAL(8):: vMod !Velocity (module)
     REAL(8):: T(1:3) !Temperature
     REAL(8):: n(1:3) !Direction of injection
+    LOGICAL:: fixDirection !The injection of particles has a fix direction defined by n
     INTEGER:: nParticles !Number of particles to introduce each time step
     CLASS(speciesGeneric), POINTER:: species !Species of injection
     INTEGER:: nEdges
@@ -69,7 +68,7 @@ MODULE moduleInject
 
   CONTAINS
     !Initialize an injection of particles
-    SUBROUTINE initInject(self, i, v, n, T, flow, units, sp, physicalSurface)
+    SUBROUTINE initInject(self, i, v, n, fixDirection, T, flow, units, sp, physicalSurface)
       USE moduleMesh
       USE moduleRefParam
       USE moduleConstParam
@@ -81,6 +80,7 @@ MODULE moduleInject
       CLASS(injectGeneric), INTENT(inout):: self
       INTEGER, INTENT(in):: i
       REAL(8), INTENT(in):: v, n(1:3), T(1:3)
+      LOGICAL, INTENT(in):: fixDirection
       INTEGER, INTENT(in):: sp, physicalSurface
       REAL(8):: tauInject
       REAL(8), INTENT(in):: flow
@@ -91,7 +91,12 @@ MODULE moduleInject
 
       self%id         =  i
       self%vMod       =  v / v_ref
-      self%n          =  n / NORM2(n)
+      IF (.NOT. fixDirection) THEN
+        self%n        =  n / NORM2(n)
+      ELSE
+        self%n        = 0.D0
+      END IF
+      self%fixDirection = fixDirection
       self%T          =  T / T_ref
       self%species    => species(sp)%obj
       tauInject = tau(self%species%n)
@@ -197,23 +202,22 @@ MODULE moduleInject
 
     END SUBROUTINE doInjects
 
-    SUBROUTINE initVelDistMaxwellian(velDist, v, T, m)
+    SUBROUTINE initVelDistMaxwellian(velDist, T, m)
       IMPLICIT NONE
 
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: v, T, m
+      REAL(8), INTENT(in):: T, m
 
-      velDist = velDistMaxwellian(v = v, vTh = DSQRT(T/m))
+      velDist = velDistMaxwellian(vTh = DSQRT(T/m))
 
     END SUBROUTINE initVelDistMaxwellian
 
-    SUBROUTINE initVelDistDelta(velDist, v)
+    SUBROUTINE initVelDistDelta(velDist)
       IMPLICIT NONE
 
       CLASS(velDistGeneric), ALLOCATABLE, INTENT(out):: velDist
-      REAL(8), INTENT(in):: v
 
-      velDist = velDistDelta(v = v)
+      velDist = velDistDelta()
 
     END SUBROUTINE initVelDistDelta
 
@@ -226,7 +230,7 @@ MODULE moduleInject
       REAL(8):: v
       v = 0.D0
 
-      v = self%v + self%vTh*randomMaxwellian()
+      v = self%vTh*randomMaxwellian()
 
     END FUNCTION randomVelMaxwellian
 
@@ -237,7 +241,7 @@ MODULE moduleInject
       CLASS(velDistDelta), INTENT(in):: self
       REAL(8):: v
 
-      v = self%v
+      v = 0.D0
 
     END FUNCTION randomVelDelta
 
@@ -256,6 +260,7 @@ MODULE moduleInject
       INTEGER:: i
       INTEGER:: n, sp
       CLASS(meshEdge), POINTER:: randomEdge
+      REAL(8):: direction(1:3)
 
       !Insert particles
       !$OMP SINGLE
@@ -300,9 +305,20 @@ MODULE moduleInject
         !Assign particle type
         partInj(n)%species => self%species
 
-        partInj(n)%v = (/ self%v(1)%obj%randomVel(), &
+        IF (self%fixDirection) THEN
-                          self%v(2)%obj%randomVel(), &
+          direction = self%n
-                          self%v(3)%obj%randomVel() /)
+
+        ELSE
+          direction = randomEdge%normal
+
+        END IF
+
+        partInj(n)%v = self%vMod*direction +  (/ self%v(1)%obj%randomVel(), &
+                                                 self%v(2)%obj%randomVel(), &
+                                                 self%v(3)%obj%randomVel() /)
+
+        print *, direction
+        print*, partInj(n)%v
 
         !Obtain natural coordinates of particle in cell
         partInj(n)%Xi = mesh%cells(partInj(n)%cell)%obj%phy2log(partInj(n)%r)