So, no idea...

Basically things do not work. I've added a correction to the node volume
in the axis which gives okays results but still this is not perfect. I
need to find a better way to do things.

Also, I've noticed that the density changes with the size of the cells,
which should not happen! I'vw to check this issue.

src/modules/mesh/2DCyl/moduleMesh2DCyl.f90

@@ -597,18 +597,18 @@ MODULE moduleMesh2DCyl
       !Computes total volume of the cell
       self%volume  = r*detJ*PI8 !2*pi * 4 (weight of 1 point 2D-Gaussian integral)
       !Computes volume per node
-      ! self%n1%v = self%n1%v + fPsi(1)*self%volume
-      ! self%n2%v = self%n2%v + fPsi(2)*self%volume
-      ! self%n3%v = self%n3%v + fPsi(3)*self%volume
-      ! self%n4%v = self%n4%v + fPsi(4)*self%volume
-      dZ = MAXVAL(self%z) - MIN(self%z)
-      r2 = MAXVAL(self%r)
-      r1 = MINVAL(self%r)
-
-      self%n1%v = self%n1%v + dZ/2.D0 * PI * (      r2**2 - 3.0D0*r1**2 + 2.0D0*r2*r1)*0.25D0
-      self%n2%v = self%n2%v + dZ/2.D0 * PI * (      r2**2 - 3.0D0*r1**2 + 2.0D0*r2*r1)*0.25D0
-      self%n3%v = self%n3%v + dZ/2.D0 * PI * (3.0D0*r2**2 -       r1**2 - 2.0D0*r2*r1)*0.25D0
-      self%n4%v = self%n4%v + dZ/2.D0 * PI * (3.0D0*r2**2 -       r1**2 - 2.0D0*r2*r1)*0.25D0
+      Xi = (/-5.D-1, -5.D-1, 0.D0/)
+      r = self%gatherF(Xi, 4, self%r)
+      self%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8
+      Xi = (/ 5.D-1, -5.D-1, 0.D0/)
+      r = self%gatherF(Xi, 4, self%r)
+      self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8
+      Xi = (/ 5.D-1,  5.D-1, 0.D0/)
+      r = self%gatherF(Xi, 4, self%r)
+      self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8
+      Xi = (/-5.D-1,  5.D-1, 0.D0/)
+      r = self%gatherF(Xi, 4, self%r)
+      self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8
 
     END SUBROUTINE volumeQuad
 

src/modules/mesh/inout/gmsh2/moduleMeshInputGmsh2.f90

@@ -296,6 +296,20 @@ MODULE moduleMeshInputGmsh2
 
       CLOSE(10)
 
+      ! Adjust node volume at axis
+      SELECT CASE(self%geometry)
+      CASE("Cyl")
+        DO n  = 1, self%numNodes
+          r = self%nodes(n)%obj%getCoordinates()
+          IF (r(2) == 0.D0) THEN
+            self%nodes(n)%obj%v = self%nodes(n)%obj%v * 3.0D0/2.0D0
+
+          END IF
+
+        END DO
+
+      END SELECT
+        
       !Call mesh connectivity
       CALL self%connectMesh
 

src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90

@@ -496,19 +496,19 @@ MODULE moduleMeshInputVTU
 
       END DO
 
-      ! ! Adjust node volume at axis
-      ! SELECT CASE(self%geometry)
-      ! CASE("Cyl")
-      !   DO n  = 1, numNodes
-      !     r = self%nodes(n)%obj%getCoordinates()
-      !     IF (r(2) == 0.D0) THEN
-      !       self%nodes(n)%obj%v = self%nodes(n)%obj%v * 2.0D0!2.0D0/3.0D0
-      !
-      !     END IF
-      !
-      !   END DO
-      !
-      ! END SELECT
+      ! Adjust node volume at axis
+      SELECT CASE(self%geometry)
+      CASE("Cyl")
+        DO n  = 1, numNodes
+          r = self%nodes(n)%obj%getCoordinates()
+          IF (r(2) == 0.D0) THEN
+            self%nodes(n)%obj%v = self%nodes(n)%obj%v * 3.0D0/2.0D0
+
+          END IF
+
+        END DO
+
+      END SELECT
         
       !Call mesh connectivity
       CALL self%connectMesh