From 6b78ae373894f0c72d746f03f4637d6c2ef7d006 Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Wed, 3 Jul 2024 22:04:28 +0200
Subject: [PATCH] Added an adhoc factor to correct injection

I think that the volume of the nodes is not being well calculated, maybe
we need a better volume calculation for this, using multiple points (as
it is done for K)
---
 src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 29 ++++++++++++++++++----
 1 file changed, 24 insertions(+), 5 deletions(-)

diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
index 7bae6b3..90773e5 100644
--- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
+++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
@@ -226,11 +226,11 @@ MODULE moduleMesh2DCyl
       REAL(8):: p1(1:2), p2(1:2)
 
       rnd = random()
-      IF (self%r(1) == 0.D0  .OR. &
-          self%r(2) == 0.D0) THEN
-        rnd = rnd**(1.D0/3.D0)
-
-      END IF
+      ! IF (self%r(1) == 0.D0  .OR. &
+      !     self%r(2) == 0.D0) THEN
+      !   rnd = rnd**(1.D0/3.D0)
+      !
+      ! END IF
 
       p1 = (/self%z(1), self%r(1) /)
       p2 = (/self%z(2), self%r(2) /)
@@ -578,6 +578,9 @@ MODULE moduleMesh2DCyl
       REAL(8):: detJ
       REAL(8):: fPsi(1:4)
       REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3)
+      REAL(8):: nodeR(1:3)
+      REAL(8), PARAMETER:: correction = 0.81D0 !TODO: No idea why this is needed, this needs more attention. Probably a better
+                                               !      calculation of the volume is needed.
 
       self%volume  = 0.D0
 
@@ -592,9 +595,25 @@ MODULE moduleMesh2DCyl
       self%volume  = r*detJ*PI8 !4*2*pi
       !Computes volume per node
       self%n1%v = self%n1%v + fPsi(1)*self%volume
+      nodeR = self%n1%getCoordinates()
+      if (nodeR(2) == 0.D0) then
+        self%n1%v = self%n1%v * correction
+      end if
       self%n2%v = self%n2%v + fPsi(2)*self%volume
+      nodeR = self%n2%getCoordinates()
+      if (nodeR(2) == 0.D0) then
+        self%n2%v = self%n2%v * correction
+      end if
       self%n3%v = self%n3%v + fPsi(3)*self%volume
+      nodeR = self%n3%getCoordinates()
+      if (nodeR(2) == 0.D0) then
+        self%n3%v = self%n3%v * correction
+      end if
       self%n4%v = self%n4%v + fPsi(4)*self%volume
+      nodeR = self%n4%getCoordinates()
+      if (nodeR(2) == 0.D0) then
+        self%n4%v = self%n4%v * correction
+      end if
 
     END SUBROUTINE volumeQuad