From 15d64f3e68f683d01de9c4a07cd2e5654169ab08 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Thu, 5 Jan 2023 21:22:13 +0100 Subject: [PATCH] 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. --- src/modules/init/moduleInput.f90 | 4 +- src/modules/mesh/0D/moduleMesh0D.f90 | 10 ++- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 56 +++++++------- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 78 ++++++++++---------- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 50 +++++++------ src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 58 ++++++++------- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 28 +++---- src/modules/mesh/moduleMesh.f90 | 18 +++-- 8 files changed, 159 insertions(+), 143 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 3ff9e5b..1a439cc 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -362,7 +362,7 @@ MODULE moduleInput nodes = mesh%cells(e)%obj%getNodes() nNodes = mesh%cells(e)%obj%nNodes ALLOCATE(fPsi(1:nNodes)) - fPsi = mesh%cells(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/)) + fPsi = mesh%cells(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/), nNodes) ALLOCATE(source(1:nNodes)) DO j = 1, nNodes source(j) = density(nodes(j)) @@ -380,7 +380,7 @@ MODULE moduleInput partNew%r = mesh%cells(e)%obj%randPos() partNew%xi = mesh%cells(e)%obj%phy2log(partNew%r) !Get mean velocity at particle position - fPsi = mesh%cells(e)%obj%fPsi(partNew%xi) + fPsi = mesh%cells(e)%obj%fPsi(partNew%xi, nNodes) DO j = 1, nNodes source(j) = velocity(nodes(j), 1) diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index a57f244..65b13a3 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -109,23 +109,25 @@ MODULE moduleMesh0D END FUNCTION randPos0D - PURE FUNCTION fPsi0D(self, Xi) RESULT(fPsi) + PURE FUNCTION fPsi0D(self, Xi, nNodes) RESULT(fPsi) IMPLICIT NONE CLASS(meshCell0D), 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.D0 END FUNCTION fPsi0D - PURE FUNCTION dPsi0D(self, Xi) RESULT(dPsi) + PURE FUNCTION dPsi0D(self, Xi, nNodes) RESULT(dPsi) IMPLICIT NONE CLASS(meshCell0D), 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 diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 0f46ba5..982c703 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -230,7 +230,7 @@ MODULE moduleMesh1DCart Xi(1) = random(-1.D0, 1.D0) Xi(2:3) = 0.D0 - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 2) r(1) = DOT_PRODUCT(fPsi, self%x) END FUNCTION randPos1DCartSegm @@ -249,7 +249,7 @@ MODULE moduleMesh1DCart self%arNodes = 0.D0 !1 point Gauss integral Xi = 0.D0 - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 2) detJ = self%detJac(Xi) l = 2.D0*detJ self%volume = l @@ -258,27 +258,29 @@ MODULE moduleMesh1DCart END SUBROUTINE areaSegm !Computes element functions at point Xi - PURE FUNCTION fPsiSegm(self, xi) RESULT(fPsi) + PURE FUNCTION fPsiSegm(self, Xi, nNodes) RESULT(fPsi) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: fPsi(1:self%nNodes) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) - fPsi(1) = 1.D0 - xi(1) - fPsi(2) = 1.D0 + xi(1) + fPsi(1) = 1.D0 - Xi(1) + fPsi(2) = 1.D0 + Xi(1) fPsi = fPsi * 5.D-1 END FUNCTION fPsiSegm !Computes element derivative shape function at Xi - PURE FUNCTION dPsiSegm(self, xi) RESULT(dPsi) + PURE FUNCTION dPsiSegm(self, Xi, nNodes) RESULT(dPsi) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: dPsi(1:3,1:self%nNodes) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) dPsi = 0.D0 @@ -314,7 +316,7 @@ MODULE moduleMesh1DCart Xi = 0.D0 DO l = 1, 3 Xi(1) = corSeg(l) - dPsi = self%dPsi(Xi) + dPsi = self%dPsi(Xi, 2) detJ = self%detJac(Xi, dPsi) invJ = self%invJac(Xi, dPsi) localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & @@ -342,7 +344,7 @@ MODULE moduleMesh1DCart DO l = 1, 3 Xi(1) = corSeg(l) detJ = self%detJac(Xi) - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 2) f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*wSeg(l)*detJ @@ -384,14 +386,14 @@ MODULE moduleMesh1DCart END FUNCTION gatherMFSegm - PURE FUNCTION insideSegm(xi) RESULT(ins) + PURE FUNCTION insideSegm(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) LOGICAL:: ins - ins = xi(1) >=-1.D0 .AND. & - xi(1) <= 1.D0 + ins = Xi(1) >=-1.D0 .AND. & + Xi(1) <= 1.D0 END FUNCTION insideSegm @@ -418,19 +420,19 @@ MODULE moduleMesh1DCart END FUNCTION phy2logSegm - !Get next element for a logical position xi - SUBROUTINE nextElementSegm(self, xi, nextElement) + !Get next element for a logical position Xi + SUBROUTINE nextElementSegm(self, Xi, nextElement) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement NULLIFY(nextElement) - IF (xi(1) < -1.D0) THEN + IF (Xi(1) < -1.D0) THEN nextElement => self%e2 - ELSEIF (xi(1) > 1.D0) THEN + ELSEIF (Xi(1) > 1.D0) THEN nextElement => self%e1 END IF @@ -440,11 +442,11 @@ MODULE moduleMesh1DCart !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS !Calculates a random position in 1D volume !Computes the element Jacobian determinant - PURE FUNCTION detJ1DCart(self, xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ1DCart(self, Xi, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshCell1DCart), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + 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 @@ -454,7 +456,7 @@ MODULE moduleMesh1DCart dPsi = dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi, 2) END IF @@ -464,11 +466,11 @@ MODULE moduleMesh1DCart END FUNCTION detJ1DCart !Computes the invers Jacobian - PURE FUNCTION invJ1DCart(self, xi, dPsi_in) RESULT(invJ) + PURE FUNCTION invJ1DCart(self, Xi, dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshCell1DCart), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + 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) @@ -478,7 +480,7 @@ MODULE moduleMesh1DCart dPsi = dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi, 2) END IF diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 8ebac17..1ce7836 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -232,7 +232,7 @@ MODULE moduleMesh1DRad Xi(1) = random(-1.D0, 1.D0) Xi(2:3) = 0.D0 - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 2) r(1) = DOT_PRODUCT(fPsi, self%r) END FUNCTION randPos1DRadSeg @@ -252,7 +252,7 @@ MODULE moduleMesh1DRad self%arNodes = 0.D0 !1 point Gauss integral Xi = 0.D0 - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 2) detJ = self%detJac(Xi) !Computes total volume of the cell r = DOT_PRODUCT(fPsi, self%r) @@ -260,38 +260,38 @@ MODULE moduleMesh1DRad self%volume = r*l !Computes volume per node Xi = (/-5.D-1, 0.D0, 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*l Xi = (/ 5.D-1, 0.D0, 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*l END SUBROUTINE areaRad !Computes element functions at point Xi - PURE FUNCTION fPsiRad(self, xi) RESULT(fPsi) + PURE FUNCTION fPsiRad(self, Xi, nNodes) RESULT(fPsi) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: fPsi(1:self%nNodes) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) - fPsi(1) = 1.D0 - xi(1) - fPsi(2) = 1.D0 + xi(1) + fPsi(1) = 1.D0 - Xi(1) + fPsi(2) = 1.D0 + Xi(1) fPsi = fPsi * 5.D-1 END FUNCTION fPsiRad !Computes element derivative shape function at Xi - PURE FUNCTION dPsiRad(self, xi) RESULT(dPsi) + PURE FUNCTION dPsiRad(self, Xi, nNodes) RESULT(dPsi) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) - REAL(8):: dPsi(1:3,1:self%nNodes) + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) dPsi = 0.D0 @@ -328,15 +328,15 @@ MODULE moduleMesh1DRad localK = 0.D0 Xi = 0.D0 DO l = 1, 3 - Xi(1) = corSeg(l) - dPsi = self%dPsi(Xi) - detJ = self%detJac(Xi, dPsi) - invJ = self%invJac(Xi, dPsi) - fPsi = self%fPsi(Xi) - r = DOT_PRODUCT(fPsi, self%r) - localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & - RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & - r*wSeg(l)/detJ + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi, 2) + detJ = self%detJac(Xi, dPsi) + invJ = self%invJac(Xi, dPsi) + fPsi = self%fPsi(Xi, 2) + r = DOT_PRODUCT(fPsi, self%r) + localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & + RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + r*wSeg(l)/detJ END DO @@ -362,7 +362,7 @@ MODULE moduleMesh1DRad DO l = 1, 3 Xi(1) = corSeg(l) detJ = self%detJac(Xi) - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 2) r = DOT_PRODUCT(fPsi, self%r) f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*r*wSeg(l)*detJ @@ -405,14 +405,14 @@ MODULE moduleMesh1DRad END FUNCTION gatherMFRad - PURE FUNCTION insideRad(xi) RESULT(ins) + PURE FUNCTION insideRad(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) LOGICAL:: ins - ins = xi(1) >=-1.D0 .AND. & - xi(1) <= 1.D0 + ins = Xi(1) >=-1.D0 .AND. & + Xi(1) <= 1.D0 END FUNCTION insideRad @@ -439,19 +439,19 @@ MODULE moduleMesh1DRad END FUNCTION phy2logRad - !Get next element for a logical position xi - SUBROUTINE nextElementRad(self, xi, nextElement) + !Get next element for a logical position Xi + SUBROUTINE nextElementRad(self, Xi, nextElement) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement NULLIFY(nextElement) - IF (xi(1) < -1.D0) THEN + IF (Xi(1) < -1.D0) THEN nextElement => self%e2 - ELSEIF (xi(1) > 1.D0) THEN + ELSEIF (Xi(1) > 1.D0) THEN nextElement => self%e1 END IF @@ -460,11 +460,11 @@ MODULE moduleMesh1DRad !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS !Computes the element Jacobian determinant - PURE FUNCTION detJ1DRad(self, xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ1DRad(self, Xi, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshCell1DRad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + 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 @@ -474,7 +474,7 @@ MODULE moduleMesh1DRad dPsi = dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi, 2) END IF @@ -484,11 +484,11 @@ MODULE moduleMesh1DRad END FUNCTION detJ1DRad !Computes the invers Jacobian - PURE FUNCTION invJ1DRad(self, xi, dPsi_in) RESULT(invJ) + PURE FUNCTION invJ1DRad(self, Xi, dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshCell1DRad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + 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):: dx(1) @@ -498,7 +498,7 @@ MODULE moduleMesh1DRad dPsi = dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi, 2) END IF diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index b3f65e4..649b8fb 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -297,19 +297,20 @@ MODULE moduleMesh2DCart !2D 1 point Gauss Quad Integral Xi = 0.D0 detJ = self%detJac(Xi)*4.D0 !4 - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 4) self%volume = detJ self%arNodes = fPsi*detJ 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(meshCell2DCartQuad), 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)) @@ -321,12 +322,13 @@ MODULE moduleMesh2DCart 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(meshCell2DCartQuad), 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 @@ -373,7 +375,7 @@ MODULE moduleMesh2DCart 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%x) r(2) = DOT_PRODUCT(fPsi, self%y) @@ -399,8 +401,8 @@ MODULE moduleMesh2DCart 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) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & @@ -431,7 +433,7 @@ MODULE moduleMesh2DCart DO m = 1, 3 Xi(2) = corQuad(m) detJ = self%detJac(Xi) - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 4) f = DOT_PRODUCT(fPsi,source) localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ @@ -521,9 +523,9 @@ MODULE moduleMesh2DCart XiO = 0.D0 DO WHILE(conv > 1.D-2) - dPsi = self%dPsi(XiO) + dPsi = self%dPsi(XiO, 4) invJ = self%invJac(XiO, dPsi) - fPsi = self%fPsi(XiO) + fPsi = self%fPsi(XiO, 4) f = (/ DOT_PRODUCT(fPsi,self%x), & DOT_PRODUCT(fPsi,self%y), & 0.D0 /) @@ -618,7 +620,7 @@ MODULE moduleMesh2DCart 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%x) r(2) = DOT_PRODUCT(fPsi, self%y) @@ -640,19 +642,20 @@ MODULE moduleMesh2DCart !2D 1 point Gauss Quad Integral Xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /) detJ = self%detJac(Xi)/2.D0 - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 4) self%volume = detJ self%arNodes = fPsi*detJ 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(meshCell2DCartTria), 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) @@ -661,12 +664,13 @@ MODULE moduleMesh2DCart 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(meshCell2DCartTria), 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 @@ -706,10 +710,10 @@ MODULE moduleMesh2DCart 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) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ END DO @@ -735,7 +739,7 @@ MODULE moduleMesh2DCart Xi(1) = Xi1Tria(l) Xi(2) = Xi2Tria(l) detJ = self%detJac(Xi) - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 4) f = DOT_PRODUCT(fPsi,source) localF = localF + f*fPsi*wTria(l)*detJ @@ -818,7 +822,7 @@ MODULE moduleMesh2DCart !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) + dPsi = self%dPsi(Xi, 4) invJ = self%invJac(Xi, dPsi) detJ = self%detJac(Xi, dPsi) Xi = MATMUL(invJ,deltaR)/detJ @@ -864,7 +868,7 @@ MODULE moduleMesh2DCart dPsi = dPsi_in ELSE - dPsi = self%dPsi(Xi) + dPsi = self%dPsi(Xi, 4) END IF @@ -889,7 +893,7 @@ MODULE moduleMesh2DCart dPsi=dPsi_in ELSE - dPsi = self%dPsi(Xi) + dPsi = self%dPsi(Xi, 4) END IF diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 8f1ee5f..adf5d4f 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -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 diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 9bd6468..3018989 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -270,7 +270,7 @@ MODULE moduleMesh3DCart !Assign proportional volume to each node Xi = (/0.25D0, 0.25D0, 0.25D0/) - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 4) volNodes = fPsi*self%volume self%n1%v = self%n1%v + volNodes(1) self%n2%v = self%n2%v + volNodes(2) @@ -298,7 +298,7 @@ MODULE moduleMesh3DCart Xi(2) = random( 0.D0, 1.D0 - Xi(1)) Xi(3) = random( 0.D0, 1.D0 - Xi(1) - Xi(2)) - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 4) r = (/ DOT_PRODUCT(fPsi, self%x), & DOT_PRODUCT(fPsi, self%y), & @@ -320,12 +320,13 @@ MODULE moduleMesh3DCart END SUBROUTINE volumeTetra !Computes element functions in point Xi - PURE FUNCTION fPsiTetra(self, Xi) RESULT(fPsi) + PURE FUNCTION fPsiTetra(self, Xi, nNodes) RESULT(fPsi) IMPLICIT NONE CLASS(meshCell3DCartTetra), 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) - Xi(3) fPsi(2) = Xi(1) @@ -335,12 +336,13 @@ MODULE moduleMesh3DCart END FUNCTION fPsiTetra !Derivative element function at coordinates Xi - PURE FUNCTION dPsiTetra(self, Xi) RESULT(dPsi) + PURE FUNCTION dPsiTetra(self, Xi, nNodes) RESULT(dPsi) IMPLICIT NONE CLASS(meshCell3DCartTetra), 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 @@ -424,10 +426,10 @@ MODULE moduleMesh3DCart Xi = 0.D0 !TODO: One point Gauss integral. Upgrade when possible Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) - 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) localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ END FUNCTION elemKTetra @@ -445,9 +447,9 @@ MODULE moduleMesh3DCart localF = 0.D0 Xi = 0.D0 Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) - dPsi = self%dPsi(Xi) + dPsi = self%dPsi(Xi, 4) detJ = self%detJac(Xi, dPsi) - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, 4) f = DOT_PRODUCT(fPsi, source) localF = f*fPsi*1.D0*detJ @@ -530,7 +532,7 @@ MODULE moduleMesh3DCart Xi = 0.D0 deltaR = (/r(1) - self%x(1), r(2) - self%y(1), r(3) - self%z(1) /) - dPsi = self%dPsi(Xi) + dPsi = self%dPsi(Xi, 4) invJ = self%invJac(Xi, dPsi) detJ = self%detJac(Xi, dPsi) Xi = MATMUL(invJ, deltaR)/detJ @@ -579,7 +581,7 @@ MODULE moduleMesh3DCart dPsi = dPsi_in ELSE - dPsi = self%dPsi(Xi) + dPsi = self%dPsi(Xi, 4) END IF @@ -604,7 +606,7 @@ MODULE moduleMesh3DCart dPsi=dPsi_in ELSE - dPsi = self%dPsi(Xi) + dPsi = self%dPsi(Xi, 4) END IF diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index d0a03f4..aeba4ba 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -215,19 +215,21 @@ MODULE moduleMesh END FUNCTION getNodesVol_interface - PURE FUNCTION fPsi_interface(self, Xi) RESULT(fPsi) + PURE FUNCTION fPsi_interface(self, Xi, nNodes) RESULT(fPsi) IMPORT:: meshCell CLASS(meshCell), 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) END FUNCTION fPsi_interface - PURE FUNCTION dPsi_interface(self, Xi) RESULT(dPsi) + PURE FUNCTION dPsi_interface(self, Xi, nNodes) RESULT(dPsi) IMPORT:: meshCell CLASS(meshCell), 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) END FUNCTION dPsi_interface @@ -530,7 +532,7 @@ MODULE moduleMesh REAL(8):: f REAL(8):: fPsi(1:self%nNodes) - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, self%nNodes) f = DOT_PRODUCT(fPsi, valNodes) END FUNCTION gatherF_scalar @@ -546,7 +548,7 @@ MODULE moduleMesh REAL(8):: f(1:n) REAL(8):: fPsi(1:self%nNodes) - fPsi = self%fPsi(Xi) + fPsi = self%fPsi(Xi, self%nNodes) f = MATMUL(fPsi, valNodes) END FUNCTION gatherF_array @@ -563,7 +565,7 @@ MODULE moduleMesh REAL(8):: dPsiR(1:3, 1:self%nNodes) REAL(8):: invJ(1:3, 1:3), detJ - dPsi = self%dPsi(Xi) + dPsi = self%dPsi(Xi, self%nNodes) detJ = self%detJac(Xi, dPsi) invJ = self%invJac(Xi, dPsi) dPsiR = MATMUL(invJ, dPsi)/detJ @@ -590,7 +592,7 @@ MODULE moduleMesh CLASS(meshNode), POINTER:: node cellNodes = self%getNodes() - fPsi = self%fPsi(part%Xi) + fPsi = self%fPsi(part%Xi, self%nNodes) tensorS = outerProduct(part%v, part%v)