From e23fc2fc2cf7d992cbc7af276bde82977462f400 Mon Sep 17 00:00:00 2001
From: Jorge Gonzalez <jorge.gonzalez@upm.es>
Date: Thu, 11 Jul 2024 11:21:38 +0200
Subject: [PATCH] Small progress

I made some small changes to how things are calculated.

I have also discovered that the issue with different density when
changing injection is not related with the node volume but with the way
injection is carried out. When loading particles from a file, all
provide the same density regardless the cell (node) volume.

I am doing testing in 2DCart as it is easier to set up.
---
 src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 | 2 +-
 src/modules/mesh/moduleMesh.f90                    | 8 +++++---
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
index 6286cfc..8637b7e 100644
--- a/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90
@@ -209,7 +209,7 @@ MODULE moduleMeshOutputVTU
       WRITE(fileID,"(8X,A)") '<CellData>'
       !Electric field
       WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="Electric Field (V m^-1)" NumberOfComponents="3">'
-      WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
+      WRITE(fileID,"(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
       WRITE(fileID,"(10X,A)") '</DataArray>'
       WRITE(fileID,"(8X,A)") '</CellData>'
 
diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90
index a2eaa5f..7390bca 100644
--- a/src/modules/mesh/moduleMesh.f90
+++ b/src/modules/mesh/moduleMesh.f90
@@ -613,6 +613,7 @@ MODULE moduleMesh
       INTEGER:: sp
       INTEGER:: i
       CLASS(meshNode), POINTER:: node
+      REAL(8):: pFraction !Particle fraction
 
       cellNodes = self%getNodes(nNodes)
       fPsi = self%fPsi(part%Xi, nNodes)
@@ -623,10 +624,11 @@ MODULE moduleMesh
 
       DO i = 1, nNodes
         node => mesh%nodes(cellNodes(i))%obj
+        pFraction = fPsi(i)*part%weight
         CALL OMP_SET_LOCK(node%lock)
-        node%output(sp)%den          = node%output(sp)%den          + part%weight*fPsi(i)
-        node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + part%weight*fPsi(i)*part%v(:)
-        node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + part%weight*fPsi(i)*tensorS
+        node%output(sp)%den          = node%output(sp)%den          + pFraction
+        node%output(sp)%mom(:)       = node%output(sp)%mom(:)       + pFraction*part%v(:)
+        node%output(sp)%tensorS(:,:) = node%output(sp)%tensorS(:,:) + pFraction*tensorS
         CALL OMP_UNSET_LOCK(node%lock)
 
       END DO