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

Small improvement for 2DCyl

Browse source 
Nothing important, but overhead in dPsi has been reduced.
This commit is contained in:
Jorge Gonzalez 2023-01-05 18:47:33 +01:00
commit 26bd73597d
2 changed files with 153 additions and 249 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -65,14 +65,13 @@ MODULE moduleMesh2DCart
!Connectivity to adjacent elements
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
REAL(8):: arNodes(1:4) = 0.D0
CONTAINS
PROCEDURE, PASS:: init => initCellQuad2DCart
PROCEDURE, PASS:: randPos => randPosCellQuad
PROCEDURE, PASS:: area => areaQuad
PROCEDURE, PASS:: fPsi => fPsiQuad
PROCEDURE, PASS:: dPsi => dPsiQuad
PROCEDURE, NOPASS, PRIVATE:: dPsiXi1 => dPsiQuadXi1
PROCEDURE, NOPASS, PRIVATE:: dPsiXi2 => dPsiQuadXi2
PROCEDURE, PASS:: fPsi => fPsiQuad
PROCEDURE, PASS:: dPsi => dPsiQuad
PROCEDURE, PASS, PRIVATE:: partialDer => partialDerQuad
PROCEDURE, PASS:: elemK => elemKQuad
PROCEDURE, PASS:: elemF => elemFQuad
@ -99,11 +98,9 @@ MODULE moduleMesh2DCart
PROCEDURE, PASS:: init => initCellTria2DCart
PROCEDURE, PASS:: randPos => randPosCellTria
PROCEDURE, PASS:: area => areaTria
PROCEDURE, PASS:: fPsi => fPsiTria
PROCEDURE, PASS:: dPsi => dPsiTria
PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1
PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2
PROCEDURE, PASS:: partialDer => partialDerTria
PROCEDURE, PASS:: fPsi => fPsiTria
PROCEDURE, PASS:: dPsi => dPsiTria
PROCEDURE, PASS, PRIVATE:: partialDer => partialDerTria
PROCEDURE, PASS:: elemK => elemKTria
PROCEDURE, PASS:: elemF => elemFTria
PROCEDURE, PASS:: gatherElectricField => gatherEFTria
@ -196,28 +193,6 @@ MODULE moduleMesh2DCart
END SUBROUTINE initEdge2DCart
!Random position in quadrilateral volume
FUNCTION randPosCellQuad(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4)
Xi(1) = random(-1.D0, 1.D0)
Xi(2) = random(-1.D0, 1.D0)
Xi(3) = 0.D0
fPsi = self%fPsi(Xi)
r(1) = DOT_PRODUCT(fPsi, self%x)
r(2) = DOT_PRODUCT(fPsi, self%y)
r(3) = 0.D0
END FUNCTION randposCellQuad
!Get nodes from edge
PURE FUNCTION getNodes2DCart(self) RESULT(n)
IMPLICIT NONE
@ -225,6 +200,7 @@ MODULE moduleMesh2DCart
CLASS(meshEdge2DCart), INTENT(in):: self
INTEGER, ALLOCATABLE:: n(:)
ALLOCATE(n(1:2))
n = (/self%n1%n, self%n2%n /)
END FUNCTION getNodes2DCart
@ -354,41 +330,21 @@ MODULE moduleMesh2DCart
dPsi = 0.D0
dPsi(1,:) = dPsiQuadXi1(Xi(2))
dPsi(2,:) = dPsiQuadXi2(Xi(1))
dPsi(1,:) = (/ -(1.D0 - Xi(2)), &
(1.D0 - Xi(2)), &
(1.D0 + Xi(2)), &
-(1.D0 + Xi(2)) /)
dPsi(2,:) = (/ -(1.D0 - Xi(1)), &
-(1.D0 + Xi(1)), &
(1.D0 + Xi(1)), &
(1.D0 - Xi(1)) /)
dPsi = dPsi * 0.25D0
END FUNCTION dPsiQuad
!Derivative element function (Xi1)
PURE FUNCTION dPsiQuadXi1(Xi2) RESULT(dPsiXi1)
IMPLICIT NONE
REAL(8),INTENT(in):: Xi2
REAL(8):: dPsiXi1(1:4)
dPsiXi1(1) = -(1.D0-Xi2)
dPsiXi1(2) = (1.D0-Xi2)
dPsiXi1(3) = (1.D0+Xi2)
dPsiXi1(4) = -(1.D0+Xi2)
dPsiXi1 = dPsiXi1*0.25D0
END FUNCTION dPsiQuadXi1
!Derivative element function (Xi2)
PURE FUNCTION dPsiQuadXi2(Xi1) RESULT(dPsiXi2)
IMPLICIT NONE
REAL(8),INTENT(in):: Xi1
REAL(8):: dPsiXi2(1:4)
dPsiXi2(1) = -(1.D0-Xi1)
dPsiXi2(2) = -(1.D0+Xi1)
dPsiXi2(3) = (1.D0+Xi1)
dPsiXi2(4) = (1.D0-Xi1)
dPsiXi2 = dPsiXi2*0.25D0
END FUNCTION dPsiQuadXi2
!Partial derivative in global coordinates
PURE SUBROUTINE partialDerQuad(self, dPsi, dx, dy)
IMPLICIT NONE
@ -396,13 +352,35 @@ MODULE moduleMesh2DCart
REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes)
REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy
dx(1) = DOT_PRODUCT(dPsi(1,:),self%x)
dx(2) = DOT_PRODUCT(dPsi(2,:),self%x)
dy(1) = DOT_PRODUCT(dPsi(1,:),self%y)
dy(2) = DOT_PRODUCT(dPsi(2,:),self%y)
dx = (/ DOT_PRODUCT(dPsi(1,1:4),self%x(1:4)), &
DOT_PRODUCT(dPsi(2,1:4),self%x(1:4)) /)
dy = (/ DOT_PRODUCT(dPsi(1,1:4),self%y(1:4)), &
DOT_PRODUCT(dPsi(2,1:4),self%y(1:4)) /)
END SUBROUTINE partialDerQuad
!Random position in quadrilateral volume
FUNCTION randPosCellQuad(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
CLASS(meshCell2DCartQuad), INTENT(in):: self
REAL(8):: r(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4)
Xi(1) = random(-1.D0, 1.D0)
Xi(2) = random(-1.D0, 1.D0)
Xi(3) = 0.D0
fPsi = self%fPsi(Xi)
r(1) = DOT_PRODUCT(fPsi, self%x)
r(2) = DOT_PRODUCT(fPsi, self%y)
r(3) = 0.D0
END FUNCTION randPosCellQuad
!Computes element local stiffness matrix
PURE FUNCTION elemKQuad(self) RESULT(localK)
IMPLICIT NONE
@ -419,14 +397,15 @@ MODULE moduleMesh2DCart
!Start 2D Gauss Quad Integral
DO l=1, 3
Xi(2) = corQuad(l)
dPsi(1,:) = self%dPsiXi1(Xi(2))
DO m = 1, 3
Xi(1) = corQuad(m)
dPsi(2,:) = self%dPsiXi2(Xi(1))
fPsi = self%fPsi(Xi)
detJ = self%detJac(Xi,dPsi)
invJ = self%invJac(Xi,dPsi)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wQuad(l)*wQuad(m)/detJ
Xi(1) = corQuad(m)
fPsi = self%fPsi(Xi)
dPsi = self%dPsi(Xi)
detJ = self%detJac(Xi,dPsi)
invJ = self%invJac(Xi,dPsi)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
MATMUL(invJ,dPsi))* &
wQuad(l)*wQuad(m)/detJ
END DO
END DO
@ -533,24 +512,25 @@ MODULE moduleMesh2DCart
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:2,1:2), f(1:2)
REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
REAL(8):: conv
!Iterative newton method to transform coordinates
conv=1.D0
XiO=0.D0
conv = 1.D0
XiO = 0.D0
DO WHILE(conv>1.D-3)
DO WHILE(conv > 1.D-3)
dPsi = self%dPsi(XiO)
invJ = self%invJac(XiO, dPsi)
fPsi = self%fPsi(XiO)
f(1) = DOT_PRODUCT(fPsi,self%x)-r(1)
f(2) = DOT_PRODUCT(fPsi,self%y)-r(2)
detJ = self%detJac(XiO,dPsi)
Xi(1:2)=XiO(1:2) - MATMUL(invJ, f)/detJ
conv=MAXVAL(DABS(Xi-XiO),1)
XiO=Xi
f = (/ DOT_PRODUCT(fPsi,self%x), &
DOT_PRODUCT(fPsi,self%y), &
0.D0 /)
f = f - r
Xi = XiO - MATMUL(invJ, f)/detJ
conv = MAXVAL(DABS(Xi-XiO),1)
XiO = Xi
END DO
@ -644,7 +624,7 @@ MODULE moduleMesh2DCart
r(2) = DOT_PRODUCT(fPsi, self%y)
r(3) = 0.D0
END FUNCTION randposCellTria
END FUNCTION randPosCellTria
!Calculates area for triangular element
PURE SUBROUTINE areaTria(self)
@ -690,37 +670,11 @@ MODULE moduleMesh2DCart
dPsi = 0.D0
dPsi(1,:) = dPsiTriaXi1(Xi(2))
dPsi(2,:) = dPsiTriaXi2(Xi(1))
dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /)
dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /)
END FUNCTION dPsiTria
!Derivative element function (Xi1)
PURE FUNCTION dPsiTriaXi1(Xi2) RESULT(dPsiXi1)
IMPLICIT NONE
REAL(8), INTENT(in):: Xi2
REAL(8):: dPsiXi1(1:3)
dPsiXi1(1) = -1.D0
dPsiXi1(2) = 1.D0
dPsiXi1(3) = 0.D0
END FUNCTION dPsiTriaXi1
!Derivative element function (Xi1)
PURE FUNCTION dPsiTriaXi2(Xi1) RESULT(dPsiXi2)
IMPLICIT NONE
REAL(8), INTENT(in):: Xi1
REAL(8):: dPsiXi2(1:3)
dPsiXi2(1) = -1.D0
dPsiXi2(2) = 0.D0
dPsiXi2(3) = 1.D0
END FUNCTION dPsiTriaXi2
PURE SUBROUTINE partialDerTria(self, dPsi, dx, dy)
IMPLICIT NONE
@ -728,10 +682,10 @@ MODULE moduleMesh2DCart
REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes)
REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy
dx(1) = DOT_PRODUCT(dPsi(1,:),self%x)
dx(2) = DOT_PRODUCT(dPsi(2,:),self%x)
dy(1) = DOT_PRODUCT(dPsi(1,:),self%y)
dy(2) = DOT_PRODUCT(dPsi(2,:),self%y)
dx = (/ DOT_PRODUCT(dPsi(1,:),self%x), &
DOT_PRODUCT(dPsi(2,:),self%x) /)
dy = (/ DOT_PRODUCT(dPsi(1,:),self%y), &
DOT_PRODUCT(dPsi(2,:),self%y) /)
END SUBROUTINE partialDerTria
@ -755,7 +709,7 @@ MODULE moduleMesh2DCart
dPsi = self%dPsi(Xi)
detJ = self%detJac(Xi,dPsi)
invJ = self%invJac(Xi,dPsi)
fPsi = self%fPsi(Xi)
fPsi = self%fPsi(Xi)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ
END DO
@ -781,7 +735,7 @@ MODULE moduleMesh2DCart
Xi(1) = Xi1Tria(l)
Xi(2) = Xi2Tria(l)
detJ = self%detJac(Xi)
fPsi = self%fPsi(Xi)
fPsi = self%fPsi(Xi)
f = DOT_PRODUCT(fPsi,source)
localF = localF + f*fPsi*wTria(l)*detJ
@ -902,8 +856,8 @@ MODULE moduleMesh2DCart
CLASS(meshCell2DCart), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes)
REAL(8):: dPsi(1:3,1:self%nNodes)
REAL(8):: dJ
REAL(8):: dPsi(1:3,1:self%nNodes)
REAL(8):: dx(1:2), dy(1:2)
IF(PRESENT(dPsi_in)) THEN
@ -915,6 +869,7 @@ MODULE moduleMesh2DCart
END IF
CALL self%partialDer(dPsi, dx, dy)
dJ = dx(1)*dy(2)-dx(2)*dy(1)
END FUNCTION detJ2DCart
@ -926,9 +881,9 @@ MODULE moduleMesh2DCart
CLASS(meshCell2DCart), INTENT(in):: self
REAL(8), INTENT(in):: Xi(1:3)
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes)
REAL(8):: invJ(1:3,1:3)
REAL(8):: dPsi(1:3,1:self%nNodes)
REAL(8):: dx(1:2), dy(1:2)
REAL(8):: invJ(1:3,1:3)
IF(PRESENT(dPsi_in)) THEN
dPsi=dPsi_in
@ -1351,7 +1306,7 @@ MODULE moduleMesh2DCart
!Revers the normal to point inside the domain
elemB%normal = - elemB%normal
END IF
END IF

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

@ -68,12 +68,10 @@ MODULE moduleMesh2DCyl
CONTAINS
PROCEDURE, PASS:: init => initCellQuad2DCyl
PROCEDURE, PASS:: randPos => randPosVolQuad
PROCEDURE, PASS:: randPos => randPosCellQuad
PROCEDURE, PASS:: area => areaQuad
PROCEDURE, PASS:: fPsi => fPsiQuad
PROCEDURE, PASS:: dPsi => dPsiQuad
PROCEDURE, NOPASS, PRIVATE:: dPsiXi1 => dPsiQuadXi1
PROCEDURE, NOPASS, PRIVATE:: dPsiXi2 => dPsiQuadXi2
PROCEDURE, PASS, PRIVATE:: partialDer => partialDerQuad
PROCEDURE, PASS:: elemK => elemKQuad
PROCEDURE, PASS:: elemF => elemFQuad
@ -98,12 +96,10 @@ MODULE moduleMesh2DCyl
CONTAINS
PROCEDURE, PASS:: init => initCellTria2DCyl
PROCEDURE, PASS:: randPos => randPosVolTria
PROCEDURE, PASS:: randPos => randPosCellTria
PROCEDURE, PASS:: area => areaTria
PROCEDURE, PASS:: fPsi => fPsiTria
PROCEDURE, PASS:: dPsi => dPsiTria
PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1
PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2
PROCEDURE, PASS, PRIVATE:: partialDer => partialDerTria
PROCEDURE, PASS:: elemK => elemKTria
PROCEDURE, PASS:: elemF => elemFTria
@ -183,6 +179,7 @@ MODULE moduleMesh2DCyl
self%z(2)-self%z(1) , &
0.D0 /)
self%normal = self%normal/NORM2(self%normal)
!Boundary index
self%boundary => boundary(bt)
ALLOCATE(self%fboundary(1:nSpecies))
@ -210,7 +207,6 @@ MODULE moduleMesh2DCyl
END FUNCTION getNodes2DCyl
PURE FUNCTION intersection2DCylEdge(self, r0) RESULT(r)
USE moduleMath
IMPLICIT NONE
CLASS(meshEdge2DCyl), INTENT(in):: self
@ -317,20 +313,16 @@ MODULE moduleMesh2DCyl
self%volume = r*detJ
!Computes volume per node
Xi = (/-5.D-1, -5.D-1, 0.D0/)
fPsi_node = self%fPsi(Xi)
r = DOT_PRODUCT(fPsi_node,self%r)
r = self%gatherF(Xi, self%r)
self%arNodes(1) = fPsi(1)*r*detJ
Xi = (/ 5.D-1, -5.D-1, 0.D0/)
fPsi_node = self%fPsi(Xi)
r = DOT_PRODUCT(fPsi_node,self%r)
r = self%gatherF(Xi, self%r)
self%arNodes(2) = fPsi(2)*r*detJ
Xi = (/ 5.D-1, 5.D-1, 0.D0/)
fPsi_node = self%fPsi(Xi)
r = DOT_PRODUCT(fPsi_node,self%r)
r = self%gatherF(Xi, self%r)
self%arNodes(3) = fPsi(3)*r*detJ
Xi = (/-5.D-1, 5.D-1, 0.D0/)
fPsi_node = self%fPsi(Xi)
r = DOT_PRODUCT(fPsi_node,self%r)
r = self%gatherF(Xi, self%r)
self%arNodes(4) = fPsi(4)*r*detJ
END SUBROUTINE areaQuad
@ -362,43 +354,20 @@ MODULE moduleMesh2DCyl
dPsi = 0.D0
dPsi(1,:) = dPsiQuadXi1(Xi(2))
dPsi(2,:) = dPsiQuadXi2(Xi(1))
dPsi(1,:) = (/ -(1.D0 - Xi(2)), &
(1.D0 - Xi(2)), &
(1.D0 + Xi(2)), &
-(1.D0 + Xi(2)) /)
dPsi(2,:) = (/ -(1.D0 - Xi(1)), &
-(1.D0 + Xi(1)), &
(1.D0 + Xi(1)), &
(1.D0 - Xi(1)) /)
dPsi = dPsi * 0.25D0
END FUNCTION dPsiQuad
!Derivative element function (Xi1)
PURE FUNCTION dPsiQuadXi1(Xi2) RESULT(dPsiXi1)
IMPLICIT NONE
REAL(8),INTENT(in):: Xi2
REAL(8):: dPsiXi1(1:4)
dPsiXi1(1) = -(1.D0 - Xi2)
dPsiXi1(2) = (1.D0 - Xi2)
dPsiXi1(3) = (1.D0 + Xi2)
dPsiXi1(4) = -(1.D0 + Xi2)
dPsiXi1 = dPsiXi1*0.25D0
END FUNCTION dPsiQuadXi1
!Derivative element function (Xi2)
PURE FUNCTION dPsiQuadXi2(Xi1) RESULT(dPsiXi2)
IMPLICIT NONE
REAL(8),INTENT(in):: Xi1
REAL(8):: dPsiXi2(1:4)
dPsiXi2(1) = -(1.D0 - Xi1)
dPsiXi2(2) = -(1.D0 + Xi1)
dPsiXi2(3) = (1.D0 + Xi1)
dPsiXi2(4) = (1.D0 - Xi1)
dPsiXi2 = dPsiXi2 * 0.25D0
END FUNCTION dPsiQuadXi2
!Partial derivative in global coordinates
PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr)
IMPLICIT NONE
@ -407,15 +376,15 @@ MODULE moduleMesh2DCyl
REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes)
REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
dz(1) = DOT_PRODUCT(dPsi(1,:),self%z)
dz(2) = DOT_PRODUCT(dPsi(2,:),self%z)
dr(1) = DOT_PRODUCT(dPsi(1,:),self%r)
dr(2) = DOT_PRODUCT(dPsi(2,:),self%r)
dz = (/ DOT_PRODUCT(dPsi(1,1:4),self%z(1:4)), &
DOT_PRODUCT(dPsi(2,1:4),self%z(1:4)) /)
dr = (/ DOT_PRODUCT(dPsi(1,1:4),self%r(1:4)), &
DOT_PRODUCT(dPsi(2,1:4),self%r(1:4)) /)
END SUBROUTINE partialDerQuad
!Random position in quadrilateral volume
FUNCTION randPosVolQuad(self) RESULT(r)
FUNCTION randPosCellQuad(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
@ -434,7 +403,7 @@ MODULE moduleMesh2DCyl
r(2) = DOT_PRODUCT(fPsi, self%r)
r(3) = 0.D0
END FUNCTION randposVolQuad
END FUNCTION randPosCellQuad
!Computes element local stiffness matrix
PURE FUNCTION elemKQuad(self) RESULT(localK)
@ -443,8 +412,9 @@ MODULE moduleMesh2DCyl
CLASS(meshCell2DCylQuad), INTENT(in):: self
REAL(8):: localK(1:self%nNodes,1:self%nNodes)
REAL(8):: r, Xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4), dPsi(1:3,1:4)
REAL(8):: r
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l, m
@ -453,17 +423,16 @@ MODULE moduleMesh2DCyl
!Start 2D Gauss Quad Integral
DO l=1, 3
Xi(2) = corQuad(l)
dPsi(1,:) = self%dPsiXi1(Xi(2))
DO m = 1, 3
Xi(1) = corQuad(m)
dPsi(2,:) = self%dPsiXi2(Xi(1))
fPsi = self%fPsi(Xi)
detJ = self%detJac(Xi,dPsi)
invJ = self%invJac(Xi,dPsi)
r = DOT_PRODUCT(fPsi,self%r)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
MATMUL(invJ,dPsi))* &
r*wQuad(l)*wQuad(m)/detJ
Xi(1) = corQuad(m)
fPsi = self%fPsi(Xi)
dPsi = self%dPsi(Xi)
detJ = self%detJac(Xi,dPsi)
invJ = self%invJac(Xi,dPsi)
r = DOT_PRODUCT(fPsi,self%r)
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), &
MATMUL(invJ,dPsi))* &
r*wQuad(l)*wQuad(m)/detJ
END DO
END DO
@ -479,8 +448,9 @@ MODULE moduleMesh2DCyl
CLASS(meshCell2DCylQuad), INTENT(in):: self
REAL(8), INTENT(in):: source(1:self%nNodes)
REAL(8):: localF(1:self%nNodes)
REAL(8):: r, Xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: fPsi(1:4)
REAL(8):: r
REAL(8):: detJ, f
INTEGER:: l, m
@ -574,23 +544,24 @@ MODULE moduleMesh2DCyl
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:2,1:2), f(1:2)
REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3)
REAL(8):: dPsi(1:3,1:4), fPsi(1:4)
REAL(8):: conv
!Iterative newton method to transform coordinates
conv=1.D0
XiO=0.D0
conv = 1.D0
XiO = 0.D0
DO WHILE(conv>1.D-3)
DO WHILE(conv > 1.D-3)
dPsi = self%dPsi(XiO)
invJ = self%invJac(XiO, dPsi)
detJ = self%detJac(XiO, dPsi)
fPsi = self%fPsi(XiO)
f = (/ DOT_PRODUCT(fPsi,self%z), &
DOT_PRODUCT(fPsi,self%r) /)
f = f - r(1:2)
Xi(1:2) = XiO(1:2) - MATMUL(invJ, f)/detJ
DOT_PRODUCT(fPsi,self%r), &
0.D0 /)
f = f - r
Xi = XiO - MATMUL(invJ, f)/detJ
conv = MAXVAL(DABS(Xi-XiO),1)
XiO = Xi
@ -667,7 +638,7 @@ MODULE moduleMesh2DCyl
END SUBROUTINE initCellTria2DCyl
!Random position in quadrilateral volume
FUNCTION randPosVolTria(self) RESULT(r)
FUNCTION randPosCellTria(self) RESULT(r)
USE moduleRandom
IMPLICIT NONE
@ -686,7 +657,7 @@ MODULE moduleMesh2DCyl
r(2) = DOT_PRODUCT(fPsi, self%r)
r(3) = 0.D0
END FUNCTION randposVolTria
END FUNCTION randPosCellTria
!Calculates area for triangular element
PURE SUBROUTINE areaTria(self)
@ -694,7 +665,8 @@ MODULE moduleMesh2DCyl
IMPLICIT NONE
CLASS(meshCell2DCylTria), INTENT(inout):: self
REAL(8):: r, Xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: r
REAL(8):: detJ
REAL(8):: fPsi(1:3)
@ -736,37 +708,11 @@ MODULE moduleMesh2DCyl
dPsi = 0.D0
dPsi(1,:) = dPsiTriaXi1(Xi(2))
dPsi(2,:) = dPsiTriaXi2(Xi(1))
dPsi(1,:) = (/ -1.D0, 1.D0, 0.D0 /)
dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /)
END FUNCTION dPsiTria
!Derivative element function (Xi1)
PURE FUNCTION dPsiTriaXi1(Xi2) RESULT(dPsiXi1)
IMPLICIT NONE
REAL(8), INTENT(in):: Xi2
REAL(8):: dPsiXi1(1:3)
dPsiXi1(1) = -1.D0
dPsiXi1(2) = 1.D0
dPsiXi1(3) = 0.D0
END FUNCTION dPsiTriaXi1
!Derivative element function (Xi1)
PURE FUNCTION dPsiTriaXi2(Xi1) RESULT(dPsiXi2)
IMPLICIT NONE
REAL(8), INTENT(in):: Xi1
REAL(8):: dPsiXi2(1:3)
dPsiXi2(1) = -1.D0
dPsiXi2(2) = 0.D0
dPsiXi2(3) = 1.D0
END FUNCTION dPsiTriaXi2
PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr)
IMPLICIT NONE
@ -774,10 +720,10 @@ MODULE moduleMesh2DCyl
REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes)
REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
dz(1) = DOT_PRODUCT(dPsi(1,:),self%z)
dz(2) = DOT_PRODUCT(dPsi(2,:),self%z)
dr(1) = DOT_PRODUCT(dPsi(1,:),self%r)
dr(2) = DOT_PRODUCT(dPsi(2,:),self%r)
dz = (/ DOT_PRODUCT(dPsi(1,:),self%z), &
DOT_PRODUCT(dPsi(2,:),self%z) /)
dr = (/ DOT_PRODUCT(dPsi(1,:),self%r), &
DOT_PRODUCT(dPsi(2,:),self%r) /)
END SUBROUTINE partialDerTria
@ -788,7 +734,8 @@ MODULE moduleMesh2DCyl
CLASS(meshCell2DCylTria), INTENT(in):: self
REAL(8):: localK(1:self%nNodes,1:self%nNodes)
REAL(8):: r, Xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: r
REAL(8):: fPsi(1:3), dPsi(1:3,1:3)
REAL(8):: invJ(1:3,1:3), detJ
INTEGER:: l
@ -820,7 +767,8 @@ MODULE moduleMesh2DCyl
REAL(8), INTENT(in):: source(1:self%nNodes)
REAL(8):: localF(1:self%nNodes)
REAL(8):: fPsi(1:3)
REAL(8):: r, Xi(1:3)
REAL(8):: Xi(1:3)
REAL(8):: r
REAL(8):: detJ, f
INTEGER:: l
@ -909,17 +857,17 @@ MODULE moduleMesh2DCyl
CLASS(meshCell2DCylTria), INTENT(in):: self
REAL(8), INTENT(in):: r(1:3)
REAL(8):: Xi(1:3)
REAL(8):: invJ(1:2,1:2), detJ
REAL(8):: deltaR(1:2)
REAL(8):: invJ(1:3,1:3), detJ
REAL(8):: deltaR(1:3)
REAL(8):: dPsi(1:3,1:3)
!Direct method to convert coordinates
Xi = 0.D0 !Irrelevant, required for input
deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /)
dPsi = self%dPsi(Xi)
invJ = self%invJac(Xi, dPsi)
detJ = self%detJac(Xi, dPsi)
Xi(1:2) = MATMUL(invJ,deltaR)/detJ
Xi = 0.D0
deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /)
dPsi = self%dPsi(Xi)
invJ = self%invJac(Xi, dPsi)
detJ = self%detJac(Xi, dPsi)
Xi = MATMUL(invJ,deltaR)/detJ
END FUNCTION phy2logTria
@ -967,6 +915,7 @@ MODULE moduleMesh2DCyl
END IF
CALL self%partialDer(dPsi, dz, dr)
dJ = dz(1)*dr(2)-dz(2)*dr(1)
END FUNCTION detJ2DCyl
@ -1006,10 +955,10 @@ MODULE moduleMesh2DCyl
INTEGER:: e, et
DO e = 1, self%numCells
!Connect Vol-Vol
!Connect Cell-Cell
DO et = 1, self%numCells
IF (e /= et) THEN
CALL connectVolVol(self%cells(e)%obj, self%cells(et)%obj)
CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj)
END IF
@ -1017,9 +966,9 @@ MODULE moduleMesh2DCyl
SELECT TYPE(self)
TYPE IS(meshParticles)
!Connect Vol-Edge
!Connect Cell-Edge
DO et = 1, self%numEdges
CALL connectVolEdge(self%cells(e)%obj, self%edges(et)%obj)
CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj)
END DO
@ -1030,7 +979,7 @@ MODULE moduleMesh2DCyl
END SUBROUTINE connectMesh2DCyl
!Selects type of elements to build connection
SUBROUTINE connectVolVol(elemA, elemB)
SUBROUTINE connectCellCell(elemA, elemB)
IMPLICIT NONE
CLASS(meshCell), INTENT(inout):: elemA
@ -1065,9 +1014,9 @@ MODULE moduleMesh2DCyl
END SELECT
END SUBROUTINE connectVolVol
END SUBROUTINE connectCellCell
SUBROUTINE connectVolEdge(elemA, elemB)
SUBROUTINE connectCellEdge(elemA, elemB)
IMPLICIT NONE
CLASS(meshCell), INTENT(inout):: elemA
@ -1088,7 +1037,7 @@ MODULE moduleMesh2DCyl
END SELECT
END SUBROUTINE connectVolEdge
END SUBROUTINE connectCellEdge
SUBROUTINE connectQuadQuad(elemA, elemB)
IMPLICIT NONE