Boundaries (EM and particles) can now be updated at every iteration

src/fpakc.f90

@@ -45,6 +45,9 @@ PROGRAM fpakc
 
   !$OMP SINGLE
   CALL verboseError("Calculating initial EM field...")
+  ! Update EM boundary models
+  call boundariesEM_update()
+
   !$OMP END SINGLE
   CALL doEMField()
   !$OMP END PARALLEL
@@ -62,6 +65,12 @@ PROGRAM fpakc
     ! Update global time step index
     timeStep = t
 
+    ! Update Particle boundary models
+    call boundariesParticle_update
+
+    ! Update EM boundary models
+    call boundariesEM_update()
+
     !Checks if a species needs to me moved in this iteration
     CALL solver%updatePushSpecies()
 

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

@@ -25,11 +25,11 @@ MODULE moduleMesh3DCart
     CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL()
     CONTAINS
       !meshEdge DEFERRED PROCEDURES
-      PROCEDURE, PASS::   init         => initEdge3DCartTria
+      PROCEDURE, PASS:: init          => initEdge3DCartTria
-      PROCEDURE, PASS::   getNodes     => getNodes3DCartTria
+      PROCEDURE, PASS:: getNodes      => getNodes3DCartTria
-      PROCEDURE, PASS::   intersection => intersection3DCartTria
+      PROCEDURE, PASS:: intersection  => intersection3DCartTria
-      PROCEDURE, PASS::   randPos      => randPosEdgeTria
+      PROCEDURE, PASS:: randPos       => randPosEdgeTria
-      procedure, pass::   center       => centerEdgeTria
+      procedure, pass:: center        => centerEdgeTria
       !PARTICULAR PROCEDURES
       PROCEDURE, NOPASS:: fPsi => fPsiTria
 

src/modules/mesh/moduleMesh.f90

@@ -257,8 +257,8 @@ MODULE moduleMesh
       PROCEDURE(detJac_interface),           DEFERRED, NOPASS:: detJac
       PROCEDURE(invJac_interface),           DEFERRED, NOPASS:: invJac
       !Procedures to get specific values in the node
-      PROCEDURE(gatherArray_interface),      DEFERRED, PASS::   gatherElectricField
+      PROCEDURE(gatherCellArray_interface),  DEFERRED, PASS::   gatherElectricField
-      PROCEDURE(gatherArray_interface),      DEFERRED, PASS::   gatherMagneticField
+      PROCEDURE(gatherCellArray_interface),  DEFERRED, PASS::   gatherMagneticField
       !Compute K and F to solve PDE on the mesh
       PROCEDURE(elemK_interface),            DEFERRED, PASS::   elemK
       PROCEDURE(elemF_interface),            DEFERRED, PASS::   elemF
@@ -390,13 +390,13 @@ MODULE moduleMesh
 
     END FUNCTION invJac_interface
 
-    PURE FUNCTION gatherArray_interface(self, Xi) RESULT(array)
+    PURE FUNCTION gatherCellArray_interface(self, Xi) RESULT(array)
       IMPORT:: meshCell
       CLASS(meshCell), INTENT(in):: self
       REAL(8), INTENT(in):: Xi(1:3)
       REAL(8):: array(1:3)
 
-    END FUNCTION gatherArray_interface
+    END FUNCTION gatherCellArray_interface
 
     PURE FUNCTION elemK_interface(self, nNodes) RESULT(localK)
       IMPORT:: meshCell
@@ -650,6 +650,7 @@ MODULE moduleMesh
   type, abstract, public:: boundaryParticleGeneric
     integer:: n
     character(:), allocatable:: name 
+    procedure(updateParticle_interface), pointer, pass:: update => null()
     contains
       procedure(applyParticle_interface), deferred, pass:: apply
 
@@ -698,6 +699,14 @@ MODULE moduleMesh
 
     end subroutine
 
+    ! Update the values of the particle boundary model
+    subroutine updateParticle_interface(self)
+      import boundaryParticleGeneric
+
+      class(boundaryParticleGeneric), intent(inout):: self
+
+    end subroutine updateParticle_interface
+
   end interface
 
   !Reflecting boundary
@@ -754,6 +763,7 @@ MODULE moduleMesh
   ! Ensures quasi-neutrality by changing the reflection coefficient
   type, public, extends(boundaryParticleGeneric):: boundaryQuasiNeutrality
     real(8):: alpha ! Reflection parameter
+    integer:: s_incident ! species index of the incident species
     type(meshEdgePointer), allocatable:: edges(:) !Array with edges
     contains
       procedure, pass:: apply => quasiNeutrality
@@ -867,6 +877,14 @@ MODULE moduleMesh
   !Array for boundaries
   type(boundaryParticleCont), allocatable, target:: boundariesParticle(:)
 
+  ! Update the particle boundary models
+  interface 
+    module subroutine boundariesParticle_update()
+
+    end subroutine boundariesParticle_update
+
+  end interface
+
   ! BOUNDARY ELECTROMAGNETIC DEFINITIONS
   ! Generic type for electromagnetic boundary conditions
   type, public, abstract:: boundaryEMGeneric
@@ -918,11 +936,11 @@ MODULE moduleMesh
 
     end subroutine applyEM_interface
 
-    ! Update the values of the boundary condition
+    ! Update the values of the EM boundary model
     subroutine updateEM_interface(self)
       import boundaryEMGeneric
 
-      class(boundaryEMGeneric), intent(in):: self
+      class(boundaryEMGeneric), intent(inout):: self
 
     end subroutine updateEM_interface
 
@@ -939,6 +957,7 @@ MODULE moduleMesh
 
   TYPE, EXTENDS(boundaryEMGeneric):: boundaryEMDirichletTime
     real(8):: potential
+    real(8):: timeFactor
     type(table1D):: temporalProfile
     contains
       ! boundaryEMGeneric DEFERRED PROCEDURES
@@ -973,6 +992,14 @@ MODULE moduleMesh
   !Information of charge and reference parameters for rho vector
   REAL(8), ALLOCATABLE:: qSpecies(:)
 
+  ! Update the EM boundary models
+  interface 
+    module subroutine boundariesEM_update()
+
+    end subroutine boundariesEM_update
+
+  end interface
+
   ! PHYSICAL SURFACES LINKING TO MESH ELEMENTS
   ! Link physical surface to edges
   type:: physicalSurface

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -89,6 +89,8 @@ submodule(moduleMesh) boundaryEM
 
         call self%temporalProfile%convert(1.D0/ti_ref, 1.D0)
 
+        self%update => updateDirichletTime
+
       end select
 
     END SUBROUTINE initDirichletTime
@@ -118,17 +120,42 @@ submodule(moduleMesh) boundaryEM
 
       CLASS(boundaryEMDirichletTime), INTENT(in):: self
       REAL(8), INTENT(inout):: vectorF(:)
-      REAL(8):: timeFactor
       integer:: n
 
-      timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin)
-
       DO n = 1, self%nNodes
-        self%nodes(n)%obj%emData%phi  = self%potential * timeFactor
+        self%nodes(n)%obj%emData%phi  = self%potential * self%timeFactor
         vectorF(self%nodes(n)%obj%n)  = self%nodes(n)%obj%emData%phi
 
       END DO
 
     END SUBROUTINE applyDirichletTime
 
+    module subroutine updateDirichletTime(self)
+      implicit none
+
+      class(boundaryEMGeneric), intent(inout):: self
+
+      select type(self)
+      type is(boundaryEMDirichletTime)
+        self%timeFactor = self%temporalProfile%get(DBLE(timeStep)*tauMin)
+
+      end select
+
+    end subroutine updateDirichletTime
+
+    module subroutine boundariesEM_update()
+      implicit none
+
+      integer:: b
+
+      do b = 1, nBoundariesEM
+        if (associated(boundariesEM(b)%obj%update)) then
+          call boundariesEM(b)%obj%update
+
+        end if
+
+      end do
+
+    end subroutine boundariesEM_update
+
 end submodule boundaryEM

src/modules/mesh/moduleMesh@boundaryParticle.f90

@@ -107,6 +107,8 @@ submodule(moduleMesh) boundaryParticle
         boundary%alpha = 0.d0
         allocate(boundary%edges(0))
 
+        boundary%update => quasiNeutrality_update
+
       end select
 
     end subroutine initQuasiNeutrality
@@ -328,24 +330,10 @@ submodule(moduleMesh) boundaryParticle
       implicit none
 
       class(boundaryQuasiNeutrality), intent(in):: self
-      class(meshEdge), intent(inout)::             edge
+      class(meshEdge), intent(inout):: edge
-      class(particle), intent(inout)::             part
+      class(particle), intent(inout):: part
-      real(8), allocatable:: density(:)
-      class(meshCell), pointer:: cell
-      real(8):: EF_dir
-      real(8):: alpha
 
-      alpha = 0.85d0
+      if (random() <= self%alpha) then
-
-      if (associated(edge%e1)) then
-        cell => edge%e1
-
-      else
-        cell => edge%e2
-
-      end if
-      
-      if (random() <= alpha) then
         call genericReflection(edge, part)
 
       else
@@ -355,14 +343,29 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine quasiNeutrality
 
+    subroutine quasiNeutrality_update(self)
+      implicit none
+
+      class(boundaryParticleGeneric), intent(inout):: self
+
+      select type(self)
+      type is(boundaryQuasiNeutrality)
+        self%alpha = 0.1d0
+
+        print *, self%alpha
+
+      end select
+
+    end subroutine quasiNeutrality_update
+
     ! Generic boundary conditions for internal use
     module subroutine genericReflection(edge, part)
       use moduleCaseParam
       use moduleSpecies
       implicit none
 
-      class(meshEdge), intent(inout)::        edge
+      class(meshEdge), intent(inout):: edge
-      class(particle), intent(inout)::        part
+      class(particle), intent(inout):: part
       !rp  = intersection between particle and edge
       !rpp = final position of particle
       !vpp = final velocity of particle
@@ -393,4 +396,19 @@ submodule(moduleMesh) boundaryParticle
 
     end subroutine genericTransparent
 
+    module subroutine boundariesParticle_update()
+      implicit none
+
+      integer:: b
+
+      do b = 1, nBoundariesParticle
+        if (associated(boundariesParticle(b)%obj%update)) then
+          call boundariesParticle(b)%obj%update
+
+        end if
+
+      end do
+
+  end subroutine boundariesParticle_update
+
 end submodule boundaryParticle