From c6470819e852aa619307ed802581ff76140eec9c Mon Sep 17 00:00:00 2001
From: JGonzalez <jorge.gonzalez@upm.es>
Date: Wed, 26 Jun 2024 15:11:01 +0200
Subject: [PATCH] Issue with calculating coordinates in quads

The third coordinate (unused) was causing some errors when it was
becomming too large.
---
 src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 10 +++++++---
 src/modules/mesh/2DCyl/moduleMesh2DCyl.f90   | 12 +++++++-----
 2 files changed, 14 insertions(+), 8 deletions(-)

diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
index 4341cb0..7bc94ce 100644
--- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90
+++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90
@@ -318,6 +318,8 @@ MODULE moduleMesh2DCart
       INTEGER, INTENT(in):: nNodes
       REAL(8):: fPsi(1:nNodes)
 
+      fPsi = 0.D0
+
       fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
                 (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
                 (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
@@ -508,15 +510,15 @@ MODULE moduleMesh2DCart
       conv = 1.D0
       XiO  = 0.D0
 
+      f(3) = 0.D0
       DO WHILE(conv > 1.D-4)
         dPsi = self%dPsi(XiO, 4)
         pDer = self%partialDer(4, dPsi)
         detJ = self%detJac(pDer)
         invJ = self%invJac(pDer)
         fPsi = self%fPsi(XiO, 4)
-        f = (/ DOT_PRODUCT(fPsi,self%x), &
-               DOT_PRODUCT(fPsi,self%y), &
-               0.D0 /) - r
+        f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), &
+                    DOT_PRODUCT(fPsi,self%y) /) - r(1:2)
         Xi = XiO - MATMUL(invJ, f)/detJ
         conv = MAXVAL(DABS(Xi-XiO),1)
         XiO = Xi
@@ -569,6 +571,7 @@ MODULE moduleMesh2DCart
       pDer = self%partialDer(4, dPsi)
       detJ = self%detJac(pDer)
       fPsi = self%fPsi(Xi, 4)
+
       !Compute total volume of the cell
       self%volume  = detJ*4.D0
       !Compute volume per node
@@ -762,6 +765,7 @@ MODULE moduleMesh2DCart
         pDer   = self%partialDer(3, dPsi)
         detJ   = self%detJac(pDer)
         invJ   = self%invJac(pDer)
+
         localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
 
       END DO
diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
index d4baedd..dfa3035 100644
--- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
+++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
@@ -326,6 +326,8 @@ MODULE moduleMesh2DCyl
       INTEGER, INTENT(in):: nNodes
       REAL(8):: fPsi(1:nNodes)
 
+      fPsi = 0.D0
+
       fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), &
                 (1.D0 + Xi(1)) * (1.D0 - Xi(2)), &
                 (1.D0 + Xi(1)) * (1.D0 + Xi(2)), &
@@ -496,7 +498,7 @@ MODULE moduleMesh2DCyl
 
     END FUNCTION elemFQuad
 
-    !Checks if Xi is inside the element
+    !Check if Xi is inside the element
     PURE FUNCTION insideQuad(Xi) RESULT(ins)
       IMPLICIT NONE
 
@@ -524,15 +526,15 @@ MODULE moduleMesh2DCyl
       conv = 1.D0
       XiO  = 0.D0
 
+      f(3) = 0.D0
       DO WHILE(conv > 1.D-4)
         dPsi = self%dPsi(XiO, 4)
         pDer = self%partialDer(4, dPsi)
         detJ = self%detJac(pDer)
         invJ = self%invJac(pDer)
         fPsi = self%fPsi(XiO, 4)
-        f = (/ DOT_PRODUCT(fPsi,self%z), &
-               DOT_PRODUCT(fPsi,self%r), &
-               0.D0 /) - r
+        f(1:2) = (/ DOT_PRODUCT(fPsi,self%z), &
+                    DOT_PRODUCT(fPsi,self%r) /) - r(1:2)
         Xi = XiO - MATMUL(invJ, f)/detJ
         conv = MAXVAL(DABS(Xi-XiO),1)
         XiO = Xi
@@ -553,7 +555,7 @@ MODULE moduleMesh2DCyl
 
       XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /)
       nextInt = MAXLOC(XiArray,1)
-      !Selects the higher value of directions and searches in that direction
+      !Select the higher value of directions and searches in that direction
       NULLIFY(neighbourElement)
       SELECT CASE (nextInt)
       CASE (1)
-- 
2.49.1