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

Issue with calculating coordinates in quads

Browse source 
The third coordinate (unused) was causing some errors when it was
becomming too large.
This commit is contained in:
Jorge Gonzalez 2024-06-26 15:11:01 +02:00
commit c6470819e8
2 changed files with 14 additions and 8 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -318,6 +318,8 @@ MODULE moduleMesh2DCart
INTEGER, INTENT(in):: nNodes INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes) REAL(8):: fPsi(1:nNodes)
fPsi = 0.D0
fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & 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)), &
(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 conv = 1.D0
XiO = 0.D0 XiO = 0.D0
f(3) = 0.D0
DO WHILE(conv > 1.D-4) DO WHILE(conv > 1.D-4)
dPsi = self%dPsi(XiO, 4) dPsi = self%dPsi(XiO, 4)
pDer = self%partialDer(4, dPsi) pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer) detJ = self%detJac(pDer)
invJ = self%invJac(pDer) invJ = self%invJac(pDer)
fPsi = self%fPsi(XiO, 4) fPsi = self%fPsi(XiO, 4)
f = (/ DOT_PRODUCT(fPsi,self%x), & f(1:2) = (/ DOT_PRODUCT(fPsi,self%x), &
DOT_PRODUCT(fPsi,self%y), & DOT_PRODUCT(fPsi,self%y) /) - r(1:2)
0.D0 /) - r
Xi = XiO - MATMUL(invJ, f)/detJ Xi = XiO - MATMUL(invJ, f)/detJ
conv = MAXVAL(DABS(Xi-XiO),1) conv = MAXVAL(DABS(Xi-XiO),1)
XiO = Xi XiO = Xi
@ -569,6 +571,7 @@ MODULE moduleMesh2DCart
pDer = self%partialDer(4, dPsi) pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer) detJ = self%detJac(pDer)
fPsi = self%fPsi(Xi, 4) fPsi = self%fPsi(Xi, 4)
!Compute total volume of the cell !Compute total volume of the cell
self%volume = detJ*4.D0 self%volume = detJ*4.D0
!Compute volume per node !Compute volume per node
@ -762,6 +765,7 @@ MODULE moduleMesh2DCart
pDer = self%partialDer(3, dPsi) pDer = self%partialDer(3, dPsi)
detJ = self%detJac(pDer) detJ = self%detJac(pDer)
invJ = self%invJac(pDer) invJ = self%invJac(pDer)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
END DO END DO

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

@ -326,6 +326,8 @@ MODULE moduleMesh2DCyl
INTEGER, INTENT(in):: nNodes INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes) REAL(8):: fPsi(1:nNodes)
fPsi = 0.D0
fPsi = (/ (1.D0 - Xi(1)) * (1.D0 - Xi(2)), & 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)), &
(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 END FUNCTION elemFQuad
!Checks if Xi is inside the element !Check if Xi is inside the element
PURE FUNCTION insideQuad(Xi) RESULT(ins) PURE FUNCTION insideQuad(Xi) RESULT(ins)
IMPLICIT NONE IMPLICIT NONE
@ -524,15 +526,15 @@ MODULE moduleMesh2DCyl
conv = 1.D0 conv = 1.D0
XiO = 0.D0 XiO = 0.D0
f(3) = 0.D0
DO WHILE(conv > 1.D-4) DO WHILE(conv > 1.D-4)
dPsi = self%dPsi(XiO, 4) dPsi = self%dPsi(XiO, 4)
pDer = self%partialDer(4, dPsi) pDer = self%partialDer(4, dPsi)
detJ = self%detJac(pDer) detJ = self%detJac(pDer)
invJ = self%invJac(pDer) invJ = self%invJac(pDer)
fPsi = self%fPsi(XiO, 4) fPsi = self%fPsi(XiO, 4)
f = (/ DOT_PRODUCT(fPsi,self%z), & f(1:2) = (/ DOT_PRODUCT(fPsi,self%z), &
DOT_PRODUCT(fPsi,self%r), & DOT_PRODUCT(fPsi,self%r) /) - r(1:2)
0.D0 /) - r
Xi = XiO - MATMUL(invJ, f)/detJ Xi = XiO - MATMUL(invJ, f)/detJ
conv = MAXVAL(DABS(Xi-XiO),1) conv = MAXVAL(DABS(Xi-XiO),1)
XiO = Xi XiO = Xi
@ -553,7 +555,7 @@ MODULE moduleMesh2DCyl
XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /) XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /)
nextInt = MAXLOC(XiArray,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) NULLIFY(neighbourElement)
SELECT CASE (nextInt) SELECT CASE (nextInt)
CASE (1) CASE (1)