JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Fixed!

Browse source 
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)
This commit is contained in:
Jorge Gonzalez 2025-08-03 22:14:19 +02:00
commit 3d7b1ce476
2 changed files with 58 additions and 56 deletions
Download patch file Download diff file Expand all files Collapse all files

58
src/modules/mesh/2DCart/moduleMesh2DCart.f90
View file

@ -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

56
src/modules/mesh/2DCyl/moduleMesh2DCyl.f90
View file

@ -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