From 3d7b1ce476a804d45fce35fbcb7a20d64c6d3186 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 3 Aug 2025 22:14:19 +0200 Subject: [PATCH] Fixed! So it seems that rectangles and triangles are now working properly. I have also checked the phy2log routine, now it seems a bit more complicated, but it is much clearer. Maybe in the future is worth rethinking to improve speed (specially for quad elements) --- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 58 ++++++++++---------- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 56 ++++++++++--------- 2 files changed, 58 insertions(+), 56 deletions(-) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index eb834e5..126d40c 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -495,34 +495,36 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad - !Transform physical coordinates to element coordinates + !Transform physical coordinates to element coordinates with a Taylor series PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv = 1.D0 - XiO = 0.D0 + conv = 1.D0 + Xi0 = 0.D0 + Xi(3) = 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(1:2) = (/ DOT_PRODUCT(fPsi,self%x), & - DOT_PRODUCT(fPsi,self%y) /) - r(1:2) - Xi = XiO - MATMUL(transpose(invJ), f)/detJ - conv = MAXVAL(DABS(Xi-XiO),1) - XiO = Xi + fPsi = self%fPsi(Xi0, 4) + x0(1) = dot_product(fPsi, self%x) + x0(2) = dot_product(fPsi, self%y) + deltaR = r(1:2) - x0 + dPsi = self%dPsi(Xi0, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ + conv = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1) + Xi0(1:2) = Xi(1:2) END DO @@ -680,8 +682,8 @@ MODULE moduleMesh2DCart dPsi = 0.D0 - dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /) - dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /) + dPsi(1,1:3) = (/ -1.D0, 1.D0, 0.D0 /) + dPsi(2,1:3) = (/ -1.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTria @@ -701,8 +703,6 @@ MODULE moduleMesh2DCart pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%y(1:3)), & DOT_PRODUCT(dPsi(2,1:3),self%y(1:3)) /) - pDer(3,3) = 1.D0 - END FUNCTION partialDerTria !Gather electric field at position Xi @@ -828,19 +828,19 @@ MODULE moduleMesh2DCart CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:3) + REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) + REAL(8):: dPsi(1:3,1:4) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - Xi = 0.D0 - deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /) - dPsi = self%dPsi(Xi, 3) - pDer = self%partialDer(3, dPsi) - invJ = transpose(self%invJac(pDer)) - detJ = self%detJac(pDer) - Xi = MATMUL(invJ,deltaR)/detJ + Xi(3) = 0.D0 + deltaR = (/ r(1) - self%x(1), r(2) - self%y(1) /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = MATMUL(pDerInv,deltaR)/detJ END FUNCTION phy2logTria diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index bb94440..1f5f33c 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -510,34 +510,36 @@ MODULE moduleMesh2DCyl END FUNCTION insideQuad - !Transform physical coordinates to element coordinates + !Transform physical coordinates to element coordinates with a Taylor series PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) + REAL(8):: Xi0(1:3), detJ, pDerInv(1:2,1:2), deltaR(1:2), x0(1:2) REAL(8):: dPsi(1:3,1:4), fPsi(1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates - conv = 1.D0 - XiO = 0.D0 + conv = 1.D0 + Xi0 = 0.D0 + Xi(3) = 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(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 + fPsi = self%fPsi(Xi0, 4) + x0(1) = dot_product(fPsi, self%z) + x0(2) = dot_product(fPsi, self%r) + deltaR = r(1:2) - x0 + dPsi = self%dPsi(Xi0, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = Xi0(1:2) + MATMUL(pDerInv, deltaR)/detJ + conv = MAXVAL(DABS(Xi(1:2)-Xi0(1:2)),1) + Xi0(1:2) = Xi(1:2) END DO @@ -858,19 +860,19 @@ MODULE moduleMesh2DCyl CLASS(meshCell2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) - REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:3) + REAL(8):: detJ, pDerInv(1:2,1:2), deltaR(1:2) + REAL(8):: dPsi(1:3,1:4) REAL(8):: pDer(1:3, 1:3) - REAL(8):: invJ(1:3, 1:3), detJ !Direct method to convert coordinates - Xi = 0.D0 - deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /) - dPsi = self%dPsi(Xi, 3) - pDer = self%partialDer(3, dPsi) - invJ = self%invJac(pDer) - detJ = self%detJac(pDer) - Xi = MATMUL(invJ,deltaR)/detJ + Xi(3) = 0.D0 + deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + pDerInv(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) + pDerInv(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) + Xi(1:2) = MATMUL(pDerInv,deltaR)/detJ END FUNCTION phy2logTria @@ -948,8 +950,8 @@ MODULE moduleMesh2DCyl invJ = 0.D0 - invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) - invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(2,1) /) + invJ(2, 1:2) = (/ -pDer(1,2), pDer(1,1) /) invJ(3, 3) = 1.D0 END FUNCTION invJ2DCyl