Passing nNodes as argument

It seems that this improves results as passing the size of the arrays as
an argument is better than getting it from self.
This commit is contained in:
Jorge Gonzalez 2023-01-05 21:22:13 +01:00
commit 15d64f3e68
8 changed files with 159 additions and 143 deletions

View file

@ -307,7 +307,7 @@ MODULE moduleMesh2DCyl
!2D 1 point Gauss Quad Integral
Xi = 0.D0
detJ = self%detJac(Xi)*PI8 !4*2*pi
fPsi = self%fPsi(Xi)
fPsi = self%fPsi(Xi, 4)
!Computes total volume of the cell
r = DOT_PRODUCT(fPsi,self%r)
self%volume = r*detJ
@ -328,29 +328,31 @@ MODULE moduleMesh2DCyl
END SUBROUTINE areaQuad
!Computes element functions in point Xi
PURE FUNCTION fPsiQuad(self, Xi) RESULT(fPsi)
PURE FUNCTION fPsiQuad(self, Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE
CLASS(meshCell2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: fPsi(1:self%nNodes)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
fPsi(1) = (1.D0-Xi(1)) * (1.D0-Xi(2))
fPsi(2) = (1.D0+Xi(1)) * (1.D0-Xi(2))
fPsi(3) = (1.D0+Xi(1)) * (1.D0+Xi(2))
fPsi(4) = (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)) /)
fPsi = fPsi*0.25D0
END FUNCTION fPsiQuad
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiQuad(self, Xi) RESULT(dPsi)
PURE FUNCTION dPsiQuad(self, Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE
CLASS(meshCell2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: dPsi(1:3,1:self%nNodes)
INTEGER, INTENT(in):: nNodes
REAL(8):: dPsi(1:3,1:nNodes)
dPsi = 0.D0
@ -397,7 +399,7 @@ MODULE moduleMesh2DCyl
Xi(2) = random(-1.D0, 1.D0)
Xi(3) = 0.D0
fPsi = self%fPsi(Xi)
fPsi = self%fPsi(Xi, 4)
r(1) = DOT_PRODUCT(fPsi, self%z)
r(2) = DOT_PRODUCT(fPsi, self%r)
@ -425,8 +427,8 @@ MODULE moduleMesh2DCyl
Xi(2) = corQuad(l)
DO m = 1, 3
Xi(1) = corQuad(m)
fPsi = self%fPsi(Xi)
dPsi = self%dPsi(Xi)
fPsi = self%fPsi(Xi, 4)
dPsi = self%dPsi(Xi, 4)
detJ = self%detJac(Xi,dPsi)
invJ = self%invJac(Xi,dPsi)
r = DOT_PRODUCT(fPsi,self%r)
@ -461,7 +463,7 @@ MODULE moduleMesh2DCyl
DO m = 1, 3
Xi(2) = corQuad(m)
detJ = self%detJac(Xi)
fPsi = self%fPsi(Xi)
fPsi = self%fPsi(Xi, 4)
r = DOT_PRODUCT(fPsi,self%r)
f = DOT_PRODUCT(fPsi,source)
localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ
@ -553,10 +555,10 @@ MODULE moduleMesh2DCyl
XiO = 0.D0
DO WHILE(conv > 1.D-2)
dPsi = self%dPsi(XiO)
dPsi = self%dPsi(XiO, 4)
invJ = self%invJac(XiO, dPsi)
detJ = self%detJac(XiO, dPsi)
fPsi = self%fPsi(XiO)
fPsi = self%fPsi(XiO, 4)
f = (/ DOT_PRODUCT(fPsi,self%z), &
DOT_PRODUCT(fPsi,self%r), &
0.D0 /)
@ -651,7 +653,7 @@ MODULE moduleMesh2DCyl
Xi(2) = random( 0.D0, 1.D0 - Xi(1))
Xi(3) = 0.D0
fPsi = self%fPsi(Xi)
fPsi = self%fPsi(Xi, 4)
r(1) = DOT_PRODUCT(fPsi, self%z)
r(2) = DOT_PRODUCT(fPsi, self%r)
@ -675,7 +677,7 @@ MODULE moduleMesh2DCyl
!2D 1 point Gauss Quad Integral
Xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
detJ = self%detJac(Xi)*PI !2PI*1/2
fPsi = self%fPsi(Xi)
fPsi = self%fPsi(Xi, 4)
!Computes total volume of the cell
r = DOT_PRODUCT(fPsi,self%r)
self%volume = r*detJ
@ -685,12 +687,13 @@ MODULE moduleMesh2DCyl
END SUBROUTINE areaTria
!Shape functions for triangular element
PURE FUNCTION fPsiTria(self, Xi) RESULT(fPsi)
PURE FUNCTION fPsiTria(self, Xi, nNodes) RESULT(fPsi)
IMPLICIT NONE
CLASS(meshCell2DCylTria), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: fPsi(1:self%nNodes)
INTEGER, INTENT(in):: nNodes
REAL(8):: fPsi(1:nNodes)
fPsi(1) = 1.D0 - Xi(1) - Xi(2)
fPsi(2) = Xi(1)
@ -699,12 +702,13 @@ MODULE moduleMesh2DCyl
END FUNCTION fPsiTria
!Derivative element function at coordinates Xi
PURE FUNCTION dPsiTria(self, Xi) RESULT(dPsi)
PURE FUNCTION dPsiTria(self, Xi, nNodes) RESULT(dPsi)
IMPLICIT NONE
CLASS(meshCell2DCylTria), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8):: dPsi(1:3,1:self%nNodes)
INTEGER, INTENT(in):: nNodes
REAL(8):: dPsi(1:3,1:nNodes)
dPsi = 0.D0
@ -746,10 +750,10 @@ MODULE moduleMesh2DCyl
DO l=1, 4
Xi(1) = Xi1Tria(l)
Xi(2) = Xi2Tria(l)
dPsi = self%dPsi(Xi)
dPsi = self%dPsi(Xi, 4)
detJ = self%detJac(Xi,dPsi)
invJ = self%invJac(Xi,dPsi)
fPsi = self%fPsi(Xi)
fPsi = self%fPsi(Xi, 4)
r = DOT_PRODUCT(fPsi,self%r)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ
@ -779,7 +783,7 @@ MODULE moduleMesh2DCyl
Xi(1) = Xi1Tria(l)
Xi(2) = Xi2Tria(l)
detJ = self%detJac(Xi)
fPsi = self%fPsi(Xi)
fPsi = self%fPsi(Xi, 4)
r = DOT_PRODUCT(fPsi,self%r)
f = DOT_PRODUCT(fPsi,source)
localF = localF + r*f*fPsi*wTria(l)*detJ
@ -864,7 +868,7 @@ MODULE moduleMesh2DCyl
!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)
dPsi = self%dPsi(Xi, 4)
invJ = self%invJac(Xi, dPsi)
detJ = self%detJac(Xi, dPsi)
Xi = MATMUL(invJ,deltaR)/detJ
@ -910,7 +914,7 @@ MODULE moduleMesh2DCyl
dPsi = dPsi_in
ELSE
dPsi = self%dPsi(Xi)
dPsi = self%dPsi(Xi, 4)
END IF
@ -935,7 +939,7 @@ MODULE moduleMesh2DCyl
dPsi=dPsi_in
ELSE
dPsi = self%dPsi(Xi)
dPsi = self%dPsi(Xi, 4)
END IF