From 26bd73597dacfb93332a93f55f08c2cbdfea8f6e Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Thu, 5 Jan 2023 18:47:33 +0100 Subject: [PATCH] Small improvement for 2DCyl Nothing important, but overhead in dPsi has been reduced. --- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 199 +++++++----------- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 203 +++++++------------ 2 files changed, 153 insertions(+), 249 deletions(-) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 13c5901..e7c6e6f 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -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 diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 9032d61..6a47027 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -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