diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90
index a20ffc2..d2fa524 100644
--- a/src/modules/mesh/moduleMesh@boundaryParticle.f90
+++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90
@@ -546,11 +546,10 @@ submodule(moduleMesh) boundaryParticle
       integer:: e, n, s
       integer, allocatable:: nodes(:)
       class(meshEdge), pointer:: edge
-      real(8), allocatable:: den_nodes(:)
       real(8), allocatable:: mom_nodes(:)
       class(meshNode), pointer:: node
-      real(8):: den_center, mom_center, vel_center
-      real(8):: vel_incident, vel_rest
+      real(8):: mom_center
+      real(8):: mom_incident, mom_rest
       real(8):: alpha
       
       select type(self)
@@ -558,52 +557,34 @@ submodule(moduleMesh) boundaryParticle
         do e = 1, size(self%edges)
           edge => self%edges(e)%obj
 
-          vel_incident = 0.d0
-          vel_rest     = 0.d0
+          mom_incident = 0.d0
+          mom_rest     = 0.d0
 
           nodes = edge%getNodes(edge%nNodes)
-          allocate(den_nodes(1:edge%nNodes))
           allocate(mom_nodes(1:edge%nNodes))
           do s = 1, nSpecies
-            den_center = 0.0d0
-            mom_center = 0.0d0
-            vel_center = 0.0d0
-            select type(sp => species(s)%obj)
-            type is(speciesCharged)
-              do n = 1, edge%nNodes
-                node => mesh%nodes(nodes(n))%obj
+            do n = 1, edge%nNodes
+              node => mesh%nodes(nodes(n))%obj
 
-                den_nodes(n) = node%output(s)%den
-                mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)
+              mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)
 
-              end do
+            end do
 
-              den_center = edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
-              if (den_center > 1.0d-10) then
-                mom_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
-                vel_center = mom_center / den_center
+            mom_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes)
 
-              else
-                vel_center = 0.0d0
-              
-              end if
-
-              if (s == self%s_incident) then
-                vel_incident = vel_center
-              else
-                vel_rest = vel_rest + vel_center
-              end if
-
-            end select
+            if (s == self%s_incident) then
+              mom_incident = mom_center
+            else
+              mom_rest = mom_rest + mom_center
+            end if
 
           end do
 
-          alpha = vel_rest - vel_incident
+          alpha = mom_rest + mom_incident
 
           ! Apply correction with a factor of 0.1 to avoid fast changes
           self%velocity_shift(edge%n) = self%velocity_shift(edge%n) + 1.0d-2 * alpha
 
-          deallocate(den_nodes)
           deallocate(mom_nodes)
 
         end do