From 086ce4f589aa05fae4b4866facf0223b5561024f Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Mon, 16 Mar 2026 12:42:49 +0100
Subject: [PATCH] Had to remove charge, it is working now and it only accounts
 for charged species

---
 .../mesh/moduleMesh@boundaryParticle.f90      | 46 +++++++++++--------
 1 file changed, 27 insertions(+), 19 deletions(-)

diff --git a/src/modules/mesh/moduleMesh@boundaryParticle.f90 b/src/modules/mesh/moduleMesh@boundaryParticle.f90
index 334ce7d..a20ffc2 100644
--- a/src/modules/mesh/moduleMesh@boundaryParticle.f90
+++ b/src/modules/mesh/moduleMesh@boundaryParticle.f90
@@ -549,7 +549,7 @@ submodule(moduleMesh) boundaryParticle
       real(8), allocatable:: den_nodes(:)
       real(8), allocatable:: mom_nodes(:)
       class(meshNode), pointer:: node
-      real(8):: den_center, vel_center
+      real(8):: den_center, mom_center, vel_center
       real(8):: vel_incident, vel_rest
       real(8):: alpha
       
@@ -565,32 +565,40 @@ submodule(moduleMesh) boundaryParticle
           allocate(den_nodes(1:edge%nNodes))
           allocate(mom_nodes(1:edge%nNodes))
           do s = 1, nSpecies
-            do n = 1, edge%nNodes
-              node => mesh%nodes(nodes(n))%obj
+            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
 
-              den_nodes(n) = node%output(s)%den
-              mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)
+                den_nodes(n) = node%output(s)%den
+                mom_nodes(n) = dot_product(node%output(s)%mom, edge%normal)
 
-            end do
+              end do
 
-            den_center = qSpecies(s)*edge%gatherF(edge%centerXi(), edge%nNodes, den_nodes)
-            if (den_center > 1.0d-10) then
-              vel_center = edge%gatherF(edge%centerXi(), edge%nNodes, mom_nodes) / den_center
+              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
 
-            else
-              vel_center = 0.0d0
-            
-            end if
+              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
+              if (s == self%s_incident) then
+                vel_incident = vel_center
+              else
+                vel_rest = vel_rest + vel_center
+              end if
+
+            end select
 
           end do
 
-          alpha = vel_rest + vel_incident
+          alpha = vel_rest - vel_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