Removed average and things seem to work

src/modules/init/moduleInput.f90

@@ -1351,8 +1351,6 @@ MODULE moduleInput
                 bound%nEdges = size(bound%edges)
                 allocate(bound%electricField(1:bound%nEdges))
                 bound%electricField = 0.0d0
-                allocate(bound%deltaElectricField(1:bound%nEdges))
-                bound%deltaElectricField = 0.0d0
 
               end if
 

src/modules/mesh/moduleMesh.f90

@@ -1118,7 +1118,6 @@ MODULE moduleMesh
     integer:: nEdges
     type(meshEdgePointer), allocatable:: edges(:) ! Edges included in the boundary
     real(8), allocatable:: electricField(:) ! Electric field normal to the edge that must be applied with a Neumann boundary condition
-    real(8), allocatable:: deltaElectricField(:) ! Accumulated change in E
     integer:: counter
     contains
       procedure, pass:: apply => applyFreeCurrent

src/modules/mesh/moduleMesh@boundaryEM.f90

@@ -378,17 +378,9 @@ submodule(moduleMesh) boundaryEM
       class(meshNode), pointer:: node
       real(8):: mom_center, den_center, charge_center
       real(8):: l_c ! Characteristic length
-      integer:: every
-      real(8):: every_inv
-
-      ! TODO: Make this an input parameter
-      every = 10
-      every_inv = 1.0d0 / real(every)
 
       select type(self)
       type is(boundaryEMFreeCurrent)
-        self%counter = self%counter + 1
-
         do e=1, self%nEdges
           mom_center    = 0.0d0
           den_center    = 0.0d0
@@ -405,6 +397,7 @@ submodule(moduleMesh) boundaryEM
                 node => mesh%nodes(nodes(n))%obj
 
                 ! Minus sign as we look at the values exiting the surface (opposite to the normal) 
+                ! Density needed to be reescaled as 'den' has no data about the node volume
                 mom_nodes(n) = - dot_product(node%output(s)%mom, edge%normal)/(node%v*Vol_ref*n_ref)
                 den_nodes(n) =   node%output(s)%den/(node%v*Vol_ref*n_ref)
 
@@ -417,10 +410,15 @@ submodule(moduleMesh) boundaryEM
             end do
 
             if (den_center > 1.0e-10) then
-              self%deltaElectricField(e) = self%deltaElectricField(e) + mom_center * charge_center / den_center * tauMin
+              l_c = mom_center / den_center * tauMin
+
+            else
+              l_c = 0.0d0
 
             end if
 
+            self%electricField(e) = self%electricField(e) - l_c * charge_center
+
             deallocate(nodes, mom_nodes, den_nodes)
 
           end associate
@@ -428,16 +426,6 @@ submodule(moduleMesh) boundaryEM
 
         end do
 
-        if (self%counter == every) then
-
-          self%electricField = self%electricField - self%deltaElectricField*every_inv
-
-          ! reset for the next iteration
-          self%counter = 0
-          self%deltaElectricField = 0.0d0
-
-        end if
-
       end select
 
     end subroutine updateFreeCurrent