From 0db76083ec1b9c49592818ed9b5c493644b51714 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Sun, 1 Jan 2023 12:12:06 +0100 Subject: [PATCH] fPsi no longer allocates memory I noticed that phy2logquad had a lot of overhead. Trying to reducing it by simplifying calls to fPsi, dPsi and such. The function for fPsi has been made so no memory is allocated and works under the assumption that the input array has the right size (1:numNodes) --- src/modules/init/moduleInput.f90 | 6 +- src/modules/mesh/0D/moduleMesh0D.f90 | 7 +- src/modules/mesh/1DCart/moduleMesh1DCart.f90 | 54 ++-- src/modules/mesh/1DRad/moduleMesh1DRad.f90 | 68 ++-- src/modules/mesh/2DCart/moduleMesh2DCart.f90 | 323 +++++++++---------- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 107 +++--- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 237 +++++++------- src/modules/mesh/moduleMesh.f90 | 15 +- 8 files changed, 408 insertions(+), 409 deletions(-) diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index ac40c3b..aaf4b08 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -361,14 +361,14 @@ MODULE moduleInput !Density at centroid of cell nodes = mesh%vols(e)%obj%getNodes() nNodes = SIZE(nodes) - fPsi = mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/)) + ALLOCATE(fPsi(1:nNodes)) + CALL mesh%vols(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/), fPsi) ALLOCATE(source(1:nNodes)) DO j = 1, nNodes source(j) = density(nodes(j)) END DO densityCen = DOT_PRODUCT(fPsi, source) - DEALLOCATE(fPsi) !Calculate number of particles nNewPart = INT(densityCen * (mesh%vols(e)%obj%volume*Vol_ref) / species(sp)%obj%weight) @@ -380,7 +380,7 @@ MODULE moduleInput partNew%r = mesh%vols(e)%obj%randPos() partNew%xi = mesh%vols(e)%obj%phy2log(partNew%r) !Get mean velocity at particle position - fPsi = mesh%vols(e)%obj%fPsi(partNew%xi) + CALL mesh%vols(e)%obj%fPsi(partNew%xi, fPsi) 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 8e8072d..0a14520 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -106,14 +106,13 @@ MODULE moduleMesh0D END FUNCTION randPos0D - PURE FUNCTION fPsi0D(xi) RESULT(fPsi) + PURE SUBROUTINE fPsi0D(xi, fPsi) REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8), INTENT(out):: fPsi(:) - ALLOCATE(fPsi(1:1)) fPsi = 1.D0 - END FUNCTION fPsi0D + END SUBROUTINE fPsi0D PURE FUNCTION gatherEF0D(self, xi) RESULT(EF) IMPLICIT NONE diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 0a846d9..3380cf0 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -230,13 +230,13 @@ MODULE moduleMesh1DCart CLASS(meshVol1DCartSegm), INTENT(in):: self REAL(8):: r(1:3) - REAL(8):: xii(1:3) + REAL(8):: Xi(1:3) REAL(8), ALLOCATABLE:: fPsi(:) - xii(1) = random(-1.D0, 1.D0) - xii(2:3) = 0.D0 + Xi(1) = random(-1.D0, 1.D0) + Xi(2:3) = 0.D0 - fPsi = self%fPsi(xii) + CALL self%fPsi(Xi, fPsi) r(1) = DOT_PRODUCT(fPsi, self%x) END FUNCTION randPos1DCartSeg @@ -249,36 +249,34 @@ MODULE moduleMesh1DCart REAL(8):: l !element length REAL(8):: fPsi(1:2) REAL(8):: detJ - REAL(8):: Xii(1:3) + REAL(8):: Xi(1:3) self%volume = 0.D0 self%arNodes = 0.D0 !1 point Gauss integral - Xii = 0.D0 - fPsi = self%fPsi(Xii) - detJ = self%detJac(Xii) + Xi = 0.D0 + CALL self%fPsi(Xi, fPsi) + detJ = self%detJac(Xi) l = 2.D0*detJ self%volume = l self%arNodes = fPsi*l END SUBROUTINE areaSegm - !Computes element functions at point xii - PURE FUNCTION fPsiSegm(xi) RESULT(fPsi) + !Computes element functions at point Xi + PURE SUBROUTINE fPsiSegm(xi, fPsi) IMPLICIT NONE REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:2)) + REAL(8), INTENT(out):: fPsi(:) fPsi(1) = 1.D0 - xi(1) fPsi(2) = 1.D0 + xi(1) fPsi = fPsi * 5.D-1 - END FUNCTION fPsiSegm + END SUBROUTINE fPsiSegm - !Computes element derivative shape function at Xii + !Computes element derivative shape function at Xi PURE FUNCTION dPsiSegm(xi) RESULT(dPsi) IMPLICIT NONE @@ -310,19 +308,19 @@ MODULE moduleMesh1DCart CLASS(meshVol1DCartSegm), INTENT(in):: self REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: Xii(1:3) + REAL(8):: Xi(1:3) REAL(8):: dPsi(1:1, 1:2) REAL(8):: invJ(1), detJ INTEGER:: l ALLOCATE(localK(1:2,1:2)) localK = 0.D0 - Xii = 0.D0 + Xi = 0.D0 DO l = 1, 3 - xii(1) = corSeg(l) - dPsi = self%dPsi(Xii) - detJ = self%detJac(Xii, dPsi) - invJ = self%invJac(Xii, dPsi) + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi, dPsi) + invJ = self%invJac(Xi, dPsi) localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & wSeg(l)/detJ @@ -339,17 +337,17 @@ MODULE moduleMesh1DCart REAL(8), ALLOCATABLE:: localF(:) REAL(8):: fPsi(1:2) REAL(8):: detJ, f - REAL(8):: Xii(1:3) + REAL(8):: Xi(1:3) INTEGER:: l ALLOCATE(localF(1:2)) localF = 0.D0 - Xii = 0.D0 + Xi = 0.D0 DO l = 1, 3 - xii(1) = corSeg(l) - detJ = self%detJac(Xii) - fPsi = self%fPsi(Xii) + Xi(1) = corSeg(l) + detJ = self%detJac(Xi) + CALL self%fPsi(Xi, fPsi) f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*wSeg(l)*detJ @@ -368,7 +366,7 @@ MODULE moduleMesh1DCart END FUNCTION insideSegm - !Gathers EF at position Xii + !Gathers EF at position Xi PURE FUNCTION gatherEFSegm(self, xi) RESULT(EF) IMPLICIT NONE @@ -407,7 +405,7 @@ MODULE moduleMesh1DCart MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), & self%n2%emData%B(3) /) - fPsi = self%fPsi(xi) + CALL self%fPsi(xi, fPsi) MF = MATMUL(fPsi, MF_Nodes) END FUNCTION gatherMFSegm diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index b8edfdd..7b09e5b 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -232,13 +232,13 @@ MODULE moduleMesh1DRad CLASS(meshVol1DRadSegm), INTENT(in):: self REAL(8):: r(1:3) - REAL(8):: xii(1:3) + REAL(8):: Xi(1:3) REAL(8), ALLOCATABLE:: fPsi(:) - xii(1) = random(-1.D0, 1.D0) - xii(2:3) = 0.D0 + Xi(1) = random(-1.D0, 1.D0) + Xi(2:3) = 0.D0 - fPsi = self%fPsi(xii) + CALL self%fPsi(Xi, fPsi) r(1) = DOT_PRODUCT(fPsi, self%r) END FUNCTION randPos1DRadSeg @@ -249,47 +249,47 @@ MODULE moduleMesh1DRad CLASS(meshVol1DRadSegm), INTENT(inout):: self REAL(8):: l !element length - REAL(8):: fPsi(1:2) + REAL(8):: fPsi(1:2), fPsi_node(1:2) REAL(8):: r REAL(8):: detJ - REAL(8):: Xii(1:3) + REAL(8):: Xi(1:3) self%volume = 0.D0 self%arNodes = 0.D0 !1 point Gauss integral - Xii = 0.D0 - fPsi = self%fPsi(Xii) - detJ = self%detJac(Xii) + Xi = 0.D0 + CALL self%fPsi(Xi, fPsi) + detJ = self%detJac(Xi) !Computes total volume of the cell r = DOT_PRODUCT(fPsi, self%r) l = 2.D0*detJ self%volume = r*l !Computes volume per node - Xii = (/-5.D-1, 0.D0, 0.D0/) - r = DOT_PRODUCT(self%fPsi(Xii),self%r) + Xi = (/-5.D-1, 0.D0, 0.D0/) + CALL self%fPsi(Xi, fPsi_node) + r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(1) = fPsi(1)*r*l - Xii = (/ 5.D-1, 0.D0, 0.D0/) - r = DOT_PRODUCT(self%fPsi(Xii),self%r) + Xi = (/ 5.D-1, 0.D0, 0.D0/) + CALL self%fPsi(Xi, fPsi_node) + r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(2) = fPsi(2)*r*l END SUBROUTINE areaRad - !Computes element functions at point xii - PURE FUNCTION fPsiRad(xi) RESULT(fPsi) + !Computes element functions at point Xi + PURE SUBROUTINE fPsiRad(xi, fPsi) IMPLICIT NONE REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:2)) + REAL(8), INTENT(out):: fPsi(:) fPsi(1) = 1.D0 - xi(1) fPsi(2) = 1.D0 + xi(1) fPsi = fPsi * 5.D-1 - END FUNCTION fPsiRad + END SUBROUTINE fPsiRad - !Computes element derivative shape function at Xii + !Computes element derivative shape function at Xi PURE FUNCTION dPsiRad(xi) RESULT(dPsi) IMPLICIT NONE @@ -322,7 +322,7 @@ MODULE moduleMesh1DRad CLASS(meshVol1DRadSegm), INTENT(in):: self REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: Xii(1:3) + REAL(8):: Xi(1:3) REAL(8):: dPsi(1:1, 1:2) REAL(8):: invJ(1), detJ REAL(8):: r, fPsi(1:2) @@ -330,13 +330,13 @@ MODULE moduleMesh1DRad ALLOCATE(localK(1:2, 1:2)) localK = 0.D0 - Xii = 0.D0 + Xi = 0.D0 DO l = 1, 3 - xii(1) = corSeg(l) - dPsi = self%dPsi(Xii) - detJ = self%detJac(Xii, dPsi) - invJ = self%invJac(Xii, dPsi) - fPsi = self%fPsi(Xii) + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi, dPsi) + invJ = self%invJac(Xi, dPsi) + CALL self%fPsi(Xi, fPsi) r = DOT_PRODUCT(fPsi, self%r) localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & @@ -357,17 +357,17 @@ MODULE moduleMesh1DRad REAL(8), ALLOCATABLE:: localF(:) REAL(8):: fPsi(1:2) REAL(8):: detJ, f, r - REAL(8):: Xii(1:3) + REAL(8):: Xi(1:3) INTEGER:: l ALLOCATE(localF(1:2)) localF = 0.D0 - Xii = 0.D0 + Xi = 0.D0 DO l = 1, 3 - xii(1) = corSeg(l) - detJ = self%detJac(Xii) - fPsi = self%fPsi(Xii) + Xi(1) = corSeg(l) + detJ = self%detJac(Xi) + CALL self%fPsi(Xi, fPsi) r = DOT_PRODUCT(fPsi, self%r) f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*r*wSeg(l)*detJ @@ -387,7 +387,7 @@ MODULE moduleMesh1DRad END FUNCTION insideRad - !Gathers EF at position Xii + !Gathers EF at position Xi PURE FUNCTION gatherEFRad(self, xi) RESULT(EF) IMPLICIT NONE @@ -426,7 +426,7 @@ MODULE moduleMesh1DRad MF_Nodes(1:2,3) = (/ self%n1%emData%B(3), & self%n2%emData%B(3) /) - fPsi = self%fPsi(xi) + CALL self%fPsi(xi, fPsi) MF = MATMUL(fPsi, MF_Nodes) END FUNCTION gatherMFRad diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index d01c0e4..c57cecc 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -11,8 +11,8 @@ MODULE moduleMesh2DCart REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) REAL(8), PARAMETER:: wQuad(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /) - REAL(8), PARAMETER:: xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /) - REAL(8), PARAMETER:: xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /) + REAL(8), PARAMETER:: Xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /) + REAL(8), PARAMETER:: Xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /) REAL(8), PARAMETER:: wTria(1:4) = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /) TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCart @@ -47,8 +47,8 @@ MODULE moduleMesh2DCart END TYPE meshVol2DCart ABSTRACT INTERFACE - PURE FUNCTION dPsi_interface(xi) RESULT(dPsi) - REAL(8), INTENT(in):: xi(1:3) + PURE FUNCTION dPsi_interface(Xi) RESULT(dPsi) + REAL(8), INTENT(in):: Xi(1:3) REAL(8), ALLOCATABLE:: dPsi(:,:) END FUNCTION dPsi_interface @@ -210,14 +210,14 @@ MODULE moduleMesh2DCart CLASS(meshVol2DCartQuad), INTENT(in):: self REAL(8):: r(1:3) - REAL(8):: xii(1:3) + REAL(8):: Xi(1:3) REAL(8), ALLOCATABLE:: fPsi(:) - xii(1) = random(-1.D0, 1.D0) - xii(2) = random(-1.D0, 1.D0) - xii(3) = 0.D0 + Xi(1) = random(-1.D0, 1.D0) + Xi(2) = random(-1.D0, 1.D0) + Xi(3) = 0.D0 - fPsi = self%fPsi(xii) + CALL self%fPsi(Xi, fPsi) r(1) = DOT_PRODUCT(fPsi, self%x) r(2) = DOT_PRODUCT(fPsi, self%y) @@ -319,78 +319,77 @@ MODULE moduleMesh2DCart IMPLICIT NONE CLASS(meshVol2DCartQuad), INTENT(inout):: self - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) REAL(8):: detJ REAL(8):: fPsi(1:4) self%volume = 0.D0 self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral - xi = 0.D0 - detJ = self%detJac(xi)*4.D0 !4 - fPsi = self%fPsi(xi) + Xi = 0.D0 + detJ = self%detJac(Xi)*4.D0 !4 + CALL self%fPsi(Xi, fPsi) self%volume = detJ self%arNodes = fPsi*detJ END SUBROUTINE areaQuad - !Computes element functions in point xi - PURE FUNCTION fPsiQuad(xi) RESULT(fPsi) + !Computes element functions in point Xi + PURE SUBROUTINE fPsiQuad(Xi, fPsi) IMPLICIT NONE - REAL(8),INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(out):: fPsi(:) - ALLOCATE(fPsi(1:4)) + 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) = (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 = fPsi*0.25D0 - END FUNCTION fPsiQuad + END SUBROUTINE fPsiQuad - !Derivative element function at coordinates xi - PURE FUNCTION dPsiQuad(xi) RESULT(dPsi) + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiQuad(Xi) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8), ALLOCATABLE:: dPsi(:,:) ALLOCATE(dPsi(1:2,1:4)) - dPsi(1,:) = dPsiQuadXi1(xi(2)) - dPsi(2,:) = dPsiQuadXi2(xi(1)) + dPsi(1,:) = dPsiQuadXi1(Xi(2)) + dPsi(2,:) = dPsiQuadXi2(Xi(1)) END FUNCTION dPsiQuad - !Derivative element function (xi1) - PURE FUNCTION dPsiQuadXi1(xi2) RESULT(dPsiXi1) + !Derivative element function (Xi1) + PURE FUNCTION dPsiQuadXi1(Xi2) RESULT(dPsiXi1) IMPLICIT NONE - REAL(8),INTENT(in):: xi2 + 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(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) + !Derivative element function (Xi2) + PURE FUNCTION dPsiQuadXi2(Xi1) RESULT(dPsiXi2) IMPLICIT NONE - REAL(8),INTENT(in):: xi1 + 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(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 @@ -415,24 +414,24 @@ MODULE moduleMesh2DCart CLASS(meshVol2DCartQuad), INTENT(in):: self REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:2,1:4) REAL(8):: invJ(1:2,1:2), detJ INTEGER:: l, m ALLOCATE(localK(1:4, 1:4)) localK=0.D0 - xi=0.D0 + Xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 3 - xi(2) = corQuad(l) - dPsi(1,:) = self%dPsiXi1(xi(2)) + 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) + Xi(1) = corQuad(m) + dPsi(2,:) = self%dPsiXi2(Xi(1)) + CALL self%fPsi(Xi, fPsi) + 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 @@ -447,20 +446,20 @@ MODULE moduleMesh2DCart CLASS(meshVol2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: source(1:) REAL(8), ALLOCATABLE:: localF(:) - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) REAL(8):: detJ, f INTEGER:: l, m ALLOCATE(localF(1:4)) localF = 0.D0 - xi = 0.D0 + Xi = 0.D0 DO l=1, 3 - xi(1) = corQuad(l) + Xi(1) = corQuad(l) DO m = 1, 3 - xi(2) = corQuad(m) - detJ = self%detJac(xi) - fPsi = self%fPsi(xi) + Xi(2) = corQuad(m) + detJ = self%detJac(Xi) + CALL self%fPsi(Xi, fPsi) f = DOT_PRODUCT(fPsi,source) localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ @@ -470,23 +469,23 @@ MODULE moduleMesh2DCart END FUNCTION elemFQuad !Checks if a particle is inside a quad element - PURE FUNCTION insideQuad(xi) RESULT(ins) + PURE FUNCTION insideQuad(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) .AND. & - (xi(2) >= -1.D0 .AND. xi(2) <= 1.D0) + ins = (Xi(1) >= -1.D0 .AND. Xi(1) <= 1.D0) .AND. & + (Xi(2) >= -1.D0 .AND. Xi(2) <= 1.D0) END FUNCTION insideQuad - !Gathers the electric field at position xi - PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF) + !Gathers the electric field at position Xi + PURE FUNCTION gatherEFQuad(self,Xi) RESULT(EF) IMPLICIT NONE CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8):: dPsi(1:2,1:4) REAL(8):: dPsiR(1:2,1:4)!Derivative of shpae functions in global coordinates REAL(8):: invJ(1:2,1:2), detJ @@ -498,9 +497,9 @@ MODULE moduleMesh2DCart self%n3%emData%phi, & self%n4%emData%phi /) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi,dPsi) + invJ = self%invJac(Xi,dPsi) dPsiR = MATMUL(invJ, dPsi)/detJ EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) @@ -508,11 +507,11 @@ MODULE moduleMesh2DCart END FUNCTION gatherEFQuad - PURE FUNCTION gatherMFQuad(self,xi) RESULT(MF) + PURE FUNCTION gatherMFQuad(self,Xi) RESULT(MF) IMPLICIT NONE CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8):: fPsi(1:4) REAL(8):: MF_Nodes(1:4,1:3) REAL(8):: MF(1:3) @@ -530,7 +529,7 @@ MODULE moduleMesh2DCart self%n3%emData%B(3), & self%n4%emData%B(3) /) - fPsi = self%fPsi(xi) + CALL self%fPsi(Xi, fPsi) MF = MATMUL(fPsi(:), MF_Nodes) END FUNCTION gatherMFQuad @@ -548,47 +547,47 @@ MODULE moduleMesh2DCart END FUNCTION getNodesQuad !Transforms physical coordinates to element coordinates - PURE FUNCTION phy2logQuad(self,r) RESULT(xN) + PURE FUNCTION phy2logQuad(self,r) RESULT(XiN) IMPLICIT NONE CLASS(meshVol2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) - REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2) + REAL(8):: XiN(1:3) + REAL(8):: XiO(1:3), detJ, invJ(1:2,1:2), f(1:2) REAL(8):: dPsi(1:2,1:4), fPsi(1:4) REAL(8):: conv !Iterative newton method to transform coordinates conv=1.D0 - xO=0.D0 + XiO=0.D0 DO WHILE(conv>1.D-4) - dPsi = self%dPsi(xO) - invJ = self%invJac(xO, dPsi) - fPsi = self%fPsi(xO) + dPsi = self%dPsi(XiO) + invJ = self%invJac(XiO, dPsi) + CALL self%fPsi(XiO, fPsi) f(1) = DOT_PRODUCT(fPsi,self%x)-r(1) f(2) = DOT_PRODUCT(fPsi,self%y)-r(2) - detJ = self%detJac(xO,dPsi) - xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ - conv=MAXVAL(DABS(xN-xO),1) - xO=xN + detJ = self%detJac(XiO,dPsi) + XiN(1:2)=XiO(1:2) - MATMUL(invJ, f)/detJ + conv=MAXVAL(DABS(XiN-XiO),1) + XiO=XiN END DO END FUNCTION phy2logQuad - !Gets the next element for a logical position xi - SUBROUTINE nextElementQuad(self, xi, nextElement) + !Gets the next element for a logical position Xi + SUBROUTINE nextElementQuad(self, Xi, nextElement) IMPLICIT NONE CLASS(meshVol2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:4) + REAL(8):: XiArray(1:4) INTEGER:: nextInt - xiArray = (/ -xi(2), xi(1), xi(2), -xi(1) /) - nextInt = MAXLOC(xiArray,1) + XiArray = (/ -Xi(2), Xi(1), Xi(2), -Xi(1) /) + nextInt = MAXLOC(XiArray,1) !Selects the higher value of directions and searches in that direction NULLIFY(nextElement) SELECT CASE (nextInt) @@ -649,14 +648,14 @@ MODULE moduleMesh2DCart CLASS(meshVol2DCartTria), INTENT(in):: self REAL(8):: r(1:3) - REAL(8):: xii(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8):: Xi(1:3) + REAL(8):: fPsi(1:3) - xii(1) = random( 0.D0, 1.D0) - xii(2) = random( 0.D0, 1.D0 - xii(1)) - xii(3) = 0.D0 + Xi(1) = random( 0.D0, 1.D0) + Xi(2) = random( 0.D0, 1.D0 - Xi(1)) + Xi(3) = 0.D0 - fPsi = self%fPsi(xii) + CALL self%fPsi(Xi, fPsi) r(1) = DOT_PRODUCT(fPsi, self%x) r(2) = DOT_PRODUCT(fPsi, self%y) @@ -669,55 +668,53 @@ MODULE moduleMesh2DCart IMPLICIT NONE CLASS(meshVol2DCartTria), INTENT(inout):: self - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) REAL(8):: detJ REAL(8):: fPsi(1:3) self%volume = 0.D0 self%arNodes = 0.D0 !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) + Xi = (/1.D0/3.D0, 1.D0/3.D0, 0.D0 /) + detJ = self%detJac(Xi)/2.D0 + CALL self%fPsi(Xi, fPsi) self%volume = detJ self%arNodes = fPsi*detJ END SUBROUTINE areaTria !Shape functions for triangular element - PURE FUNCTION fPsiTria(xi) RESULT(fPsi) + PURE SUBROUTINE fPsiTria(Xi, fPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(out):: fPsi(:) - ALLOCATE(fPsi(1:3)) + fPsi(1) = 1.D0 - Xi(1) - Xi(2) + fPsi(2) = Xi(1) + fPsi(3) = Xi(2) - fPsi(1) = 1.D0 - xi(1) - xi(2) - fPsi(2) = xi(1) - fPsi(3) = xi(2) + END SUBROUTINE fPsiTria - END FUNCTION fPsiTria - - !Derivative element function at coordinates xi - PURE FUNCTION dPsiTria(xi) RESULT(dPsi) + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiTria(Xi) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8), ALLOCATABLE:: dPsi(:,:) ALLOCATE(dPsi(1:2,1:3)) - dPsi(1,:) = dPsiTriaXi1(xi(2)) - dPsi(2,:) = dPsiTriaXi2(xi(1)) + dPsi(1,:) = dPsiTriaXi1(Xi(2)) + dPsi(2,:) = dPsiTriaXi2(Xi(1)) END FUNCTION dPsiTria - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi1(xi2) RESULT(dPsiXi1) + !Derivative element function (Xi1) + PURE FUNCTION dPsiTriaXi1(Xi2) RESULT(dPsiXi1) IMPLICIT NONE - REAL(8), INTENT(in):: xi2 + REAL(8), INTENT(in):: Xi2 REAL(8):: dPsiXi1(1:3) dPsiXi1(1) = -1.D0 @@ -726,11 +723,11 @@ MODULE moduleMesh2DCart END FUNCTION dPsiTriaXi1 - !Derivative element function (xi1) - PURE FUNCTION dPsiTriaXi2(xi1) RESULT(dPsiXi2) + !Derivative element function (Xi1) + PURE FUNCTION dPsiTriaXi2(Xi1) RESULT(dPsiXi2) IMPLICIT NONE - REAL(8), INTENT(in):: xi1 + REAL(8), INTENT(in):: Xi1 REAL(8):: dPsiXi2(1:3) dPsiXi2(1) = -1.D0 @@ -759,22 +756,22 @@ MODULE moduleMesh2DCart CLASS(meshVol2DCartTria), INTENT(in):: self REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:3), dPsi(1:2,1:3) REAL(8):: invJ(1:2,1:2), detJ INTEGER:: l ALLOCATE(localK(1:4, 1:4)) localK=0.D0 - xi=0.D0 + Xi=0.D0 !Start 2D Gauss Quad Integral DO l=1, 4 - xi(1) = xi1Tria(l) - xi(2) = xi2Tria(l) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) - fPsi = self%fPsi(xi) + Xi(1) = Xi1Tria(l) + Xi(2) = Xi2Tria(l) + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi,dPsi) + invJ = self%invJac(Xi,dPsi) + CALL self%fPsi(Xi, fPsi) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ END DO @@ -789,19 +786,19 @@ MODULE moduleMesh2DCart REAL(8), INTENT(in):: source(1:) REAL(8), ALLOCATABLE:: localF(:) REAL(8):: fPsi(1:3) - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) REAL(8):: detJ, f INTEGER:: l ALLOCATE(localF(1:3)) localF = 0.D0 - xi = 0.D0 + Xi = 0.D0 !Start 2D Gauss Quad Integral DO l=1, 4 - xi(1) = xi1Tria(l) - xi(2) = xi2Tria(l) - detJ = self%detJac(xi) - fPsi = self%fPsi(xi) + Xi(1) = Xi1Tria(l) + Xi(2) = Xi2Tria(l) + detJ = self%detJac(Xi) + CALL self%fPsi(Xi, fPsi) f = DOT_PRODUCT(fPsi,source) localF = localF + f*fPsi*wTria(l)*detJ @@ -809,24 +806,24 @@ MODULE moduleMesh2DCart END FUNCTION elemFTria - PURE FUNCTION insideTria(xi) RESULT(ins) + PURE FUNCTION insideTria(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) LOGICAL:: ins - ins = xi(1) >= 0.D0 .AND. & - xi(2) >= 0.D0 .AND. & - 1.D0 - xi(1) - xi(2) >= 0.D0 + ins = Xi(1) >= 0.D0 .AND. & + Xi(2) >= 0.D0 .AND. & + 1.D0 - Xi(1) - Xi(2) >= 0.D0 END FUNCTION insideTria - !Gathers the electric field at position xi - PURE FUNCTION gatherEFTria(self,xi) RESULT(EF) + !Gathers the electric field at position Xi + PURE FUNCTION gatherEFTria(self,Xi) RESULT(EF) IMPLICIT NONE CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8):: dPsi(1:2,1:3) REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates REAL(8):: invJ(1:2,1:2), detJ @@ -837,9 +834,9 @@ MODULE moduleMesh2DCart self%n2%emData%phi, & self%n3%emData%phi /) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi,dPsi) - invJ = self%invJac(xi,dPsi) + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi,dPsi) + invJ = self%invJac(Xi,dPsi) dPsiR = MATMUL(invJ, dPsi)/detJ EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) @@ -847,11 +844,11 @@ MODULE moduleMesh2DCart END FUNCTION gatherEFTria - PURE FUNCTION gatherMFTria(self,xi) RESULT(MF) + PURE FUNCTION gatherMFTria(self,Xi) RESULT(MF) IMPLICIT NONE CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8):: fPsi(1:3) REAL(8):: MF_Nodes(1:3,1:3) REAL(8):: MF(1:3) @@ -866,7 +863,7 @@ MODULE moduleMesh2DCart self%n2%emData%B(3), & self%n3%emData%B(3) /) - fPsi = self%fPsi(xi) + CALL self%fPsi(Xi, fPsi) MF = MATMUL(fPsi, MF_Nodes) END FUNCTION gatherMFTria @@ -884,37 +881,37 @@ MODULE moduleMesh2DCart END FUNCTION getNodesTria !Transforms physical coordinates to element coordinates - PURE FUNCTION phy2logTria(self,r) RESULT(xi) + PURE FUNCTION phy2logTria(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshVol2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) REAL(8):: invJ(1:2,1:2), detJ REAL(8):: deltaR(1:2) REAL(8):: dPsi(1:2,1:3) !Direct method to convert coordinates - xi = 0.D0 !Irrelevant, required for input + Xi = 0.D0 !Irrelevant, required for input deltaR = (/ r(1) - self%x(1), r(2) - self%y(1) /) - dPsi = self%dPsi(xi) - invJ = self%invJac(xi, dPsi) - detJ = self%detJac(xi, dPsi) - xi(1:2) = MATMUL(invJ,deltaR)/detJ + dPsi = self%dPsi(Xi) + invJ = self%invJac(Xi, dPsi) + detJ = self%detJac(Xi, dPsi) + Xi(1:2) = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria - SUBROUTINE nextElementTria(self, xi, nextElement) + SUBROUTINE nextElementTria(self, Xi, nextElement) IMPLICIT NONE CLASS(meshVol2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:3) + REAL(8):: XiArray(1:3) INTEGER:: nextInt - xiArray = (/ xi(2), 1.D0-xi(1)-xi(2), xi(1) /) - nextInt = MINLOC(xiArray,1) + XiArray = (/ Xi(2), 1.D0-Xi(1)-Xi(2), Xi(1) /) + nextInt = MINLOC(XiArray,1) NULLIFY(nextElement) SELECT CASE (nextInt) CASE (1) @@ -929,11 +926,11 @@ MODULE moduleMesh2DCart !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 2D !Computes element Jacobian determinant - PURE FUNCTION detJ2DCart(self, xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ2DCart(self, Xi, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshVol2DCart), 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:,1:) REAL(8), ALLOCATABLE:: dPsi(:,:) REAL(8):: dJ @@ -943,7 +940,7 @@ MODULE moduleMesh2DCart dPsi = dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi) END IF @@ -953,11 +950,11 @@ MODULE moduleMesh2DCart END FUNCTION detJ2DCart !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJ2DCart(self,xi,dPsi_in) RESULT(invJ) + PURE FUNCTION invJ2DCart(self,Xi,dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshVol2DCart), 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:,1:) REAL(8), ALLOCATABLE:: dPsi(:,:) REAL(8):: dx(1:2), dy(1:2) @@ -967,7 +964,7 @@ MODULE moduleMesh2DCart dPsi=dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi) END IF diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index c576c71..2869cb3 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -310,49 +310,52 @@ MODULE moduleMesh2DCyl CLASS(meshVol2DCylQuad), INTENT(inout):: self REAL(8):: r, xi(1:3) REAL(8):: detJ - REAL(8):: fPsi(1:4) + REAL(8):: fPsi(1:4), fPsi_node(1:4) self%volume = 0.D0 self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral xi = 0.D0 detJ = self%detJac(xi)*PI8 !4*2*pi - fPsi = self%fPsi(xi) + CALL self%fPsi(xi, fPsi) !Computes total volume of the cell r = DOT_PRODUCT(fPsi,self%r) self%volume = r*detJ !Computes volume per node xi = (/-5.D-1, -5.D-1, 0.D0/) - r = DOT_PRODUCT(self%fPsi(xi),self%r) + CALL self%fPsi(xi, fPsi_node) + r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(1) = fPsi(1)*r*detJ xi = (/ 5.D-1, -5.D-1, 0.D0/) - r = DOT_PRODUCT(self%fPsi(xi),self%r) + CALL self%fPsi(xi, fPsi_node) + r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(2) = fPsi(2)*r*detJ xi = (/ 5.D-1, 5.D-1, 0.D0/) - r = DOT_PRODUCT(self%fPsi(xi),self%r) + CALL self%fPsi(xi, fPsi_node) + r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(3) = fPsi(3)*r*detJ xi = (/-5.D-1, 5.D-1, 0.D0/) - r = DOT_PRODUCT(self%fPsi(xi),self%r) + CALL self%fPsi(xi, fPsi_node) + r = DOT_PRODUCT(fPsi_node,self%r) self%arNodes(4) = fPsi(4)*r*detJ END SUBROUTINE areaQuad !Computes element functions in point xi - PURE FUNCTION fPsiQuad(xi) RESULT(fPsi) + PURE SUBROUTINE fPsiQuad(xi, fPsi) IMPLICIT NONE - REAL(8),INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(out):: fPsi(:) - ALLOCATE(fPsi(1:4)) + 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) = (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 = fPsi*0.25D0 - END FUNCTION fPsiQuad + END SUBROUTINE fPsiQuad !Derivative element function at coordinates xi PURE FUNCTION dPsiQuad(xi) RESULT(dPsi) @@ -375,10 +378,11 @@ MODULE moduleMesh2DCyl 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(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 @@ -390,11 +394,12 @@ MODULE moduleMesh2DCyl 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 + 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 @@ -427,7 +432,7 @@ MODULE moduleMesh2DCyl xii(2) = random(-1.D0, 1.D0) xii(3) = 0.D0 - fPsi = self%fPsi(xii) + CALL self%fPsi(xii, fPsi) r(1) = DOT_PRODUCT(fPsi, self%z) r(2) = DOT_PRODUCT(fPsi, self%r) @@ -457,7 +462,7 @@ MODULE moduleMesh2DCyl DO m = 1, 3 xi(1) = corQuad(m) dPsi(2,:) = self%dPsiXi2(xi(1)) - fPsi = self%fPsi(xi) + CALL self%fPsi(xi, fPsi) detJ = self%detJac(xi,dPsi) invJ = self%invJac(xi,dPsi) r = DOT_PRODUCT(fPsi,self%r) @@ -492,7 +497,7 @@ MODULE moduleMesh2DCyl DO m = 1, 3 xi(2) = corQuad(m) detJ = self%detJac(xi) - fPsi = self%fPsi(xi) + CALL self%fPsi(xi, fPsi) r = DOT_PRODUCT(fPsi,self%r) f = DOT_PRODUCT(fPsi,source) localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ @@ -564,7 +569,7 @@ MODULE moduleMesh2DCyl self%n3%emData%B(3), & self%n4%emData%B(3) /) - fPsi = self%fPsi(xi) + CALL self%fPsi(xi, fPsi) MF = MATMUL(fPsi(:), MF_Nodes) END FUNCTION gatherMFQuad @@ -582,30 +587,30 @@ MODULE moduleMesh2DCyl END FUNCTION getNodesQuad !Transforms physical coordinates to element coordinates - PURE FUNCTION phy2logQuad(self,r) RESULT(xN) + PURE FUNCTION phy2logQuad(self,r) RESULT(XiN) IMPLICIT NONE CLASS(meshVol2DCylQuad), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) - REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2) + REAL(8):: XiN(1:3) + REAL(8):: XiO(1:3), detJ, invJ(1:2,1:2), f(1:2) REAL(8):: dPsi(1:2,1:4), fPsi(1:4) REAL(8):: conv !Iterative newton method to transform coordinates conv=1.D0 - xO=0.D0 + XiO=0.D0 - DO WHILE(conv>1.D-4) - dPsi = self%dPsi(xO) - invJ = self%invJac(xO, dPsi) - fPsi = self%fPsi(xO) - f(1) = DOT_PRODUCT(fPsi,self%z)-r(1) - f(2) = DOT_PRODUCT(fPsi,self%r)-r(2) - detJ = self%detJac(xO,dPsi) - xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ - conv=MAXVAL(DABS(xN-xO),1) - xO=xN + DO WHILE(conv>1.D-3) + CALL self%fPsi(XiO, fPsi) + f = (/ DOT_PRODUCT(fPsi,self%z)-r(1), & + DOT_PRODUCT(fPsi,self%r)-r(2) /) + dPsi = self%dPsi(XiO) + invJ = self%invJac(XiO, dPsi) + detJ = self%detJac(XiO,dPsi) + XiN(1:2)=XiO(1:2) - MATMUL(invJ, f)/detJ + conv=MAXVAL(DABS(XiN-XiO),1) + XiO=XiN END DO @@ -690,7 +695,7 @@ MODULE moduleMesh2DCyl xii(2) = random( 0.D0, 1.D0 - xii(1)) xii(3) = 0.D0 - fPsi = self%fPsi(xii) + CALL self%fPsi(xii, fPsi) r(1) = DOT_PRODUCT(fPsi, self%z) r(2) = DOT_PRODUCT(fPsi, self%r) @@ -713,7 +718,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) + CALL self%fPsi(xi, fPsi) !Computes total volume of the cell r = DOT_PRODUCT(fPsi,self%r) self%volume = r*detJ @@ -723,19 +728,17 @@ MODULE moduleMesh2DCyl END SUBROUTINE areaTria !Shape functions for triangular element - PURE FUNCTION fPsiTria(xi) RESULT(fPsi) + PURE SUBROUTINE fPsiTria(xi, fPsi) IMPLICIT NONE REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:3)) + REAL(8), INTENT(out):: fPsi(:) fPsi(1) = 1.D0 - xi(1) - xi(2) fPsi(2) = xi(1) fPsi(3) = xi(2) - END FUNCTION fPsiTria + END SUBROUTINE fPsiTria !Derivative element function at coordinates xi PURE FUNCTION dPsiTria(xi) RESULT(dPsi) @@ -813,7 +816,7 @@ MODULE moduleMesh2DCyl dPsi = self%dPsi(xi) detJ = self%detJac(xi,dPsi) invJ = self%invJac(xi,dPsi) - fPsi = self%fPsi(xi) + CALL self%fPsi(xi, fPsi) r = DOT_PRODUCT(fPsi,self%r) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ @@ -843,7 +846,7 @@ MODULE moduleMesh2DCyl xi(1) = xi1Tria(l) xi(2) = xi2Tria(l) detJ = self%detJac(xi) - fPsi = self%fPsi(xi) + CALL self%fPsi(xi, fPsi) r = DOT_PRODUCT(fPsi,self%r) f = DOT_PRODUCT(fPsi,source) localF = localF + r*f*fPsi*wTria(l)*detJ @@ -910,7 +913,7 @@ MODULE moduleMesh2DCyl self%n2%emData%B(3), & self%n3%emData%B(3) /) - fPsi = self%fPsi(xi) + CALL self%fPsi(xi, fPsi) MF = MATMUL(fPsi, MF_Nodes) END FUNCTION gatherMFTria diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 0add18a..4b29c16 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -41,8 +41,8 @@ MODULE moduleMesh3DCart END TYPE meshVol3DCart ABSTRACT INTERFACE - PURE FUNCTION dPsi_interface(xii) RESULT(dPsi) - REAL(8), INTENT(in):: xii(1:3) + PURE FUNCTION dPsi_interface(Xi) RESULT(dPsi) + REAL(8), INTENT(in):: Xi(1:3) REAL(8), ALLOCATABLE:: dPsi(:,:) END FUNCTION dPsi_interface @@ -71,8 +71,8 @@ MODULE moduleMesh3DCart PROCEDURE, PASS:: calcVol => volumeTetra PROCEDURE, NOPASS:: fPsi => fPsiTetra PROCEDURE, NOPASS:: dPsi => dPsiTetra - PROCEDURE, NOPASS:: dPsiXi1 => dPsiTetraXii1 - PROCEDURE, NOPASS:: dPsiXi2 => dPsiTetraXii2 + PROCEDURE, NOPASS:: dPsiXi1 => dPsiTetraXi1 + PROCEDURE, NOPASS:: dPsiXi2 => dPsiTetraXi2 PROCEDURE, PASS:: partialDer => partialDerTetra PROCEDURE, PASS:: elemK => elemKTetra PROCEDURE, PASS:: elemF => elemFTetra @@ -213,14 +213,14 @@ MODULE moduleMesh3DCart CLASS(meshEdge3DCartTria), INTENT(in):: self REAL(8):: r(1:3) - REAL(8):: xii(1:3) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:3) - xii(1) = random( 0.D0, 1.D0) - xii(2) = random( 0.D0, 1.D0 - xii(1)) - xii(3) = 0.D0 + Xi(1) = random( 0.D0, 1.D0) + Xi(2) = random( 0.D0, 1.D0 - Xi(1)) + Xi(3) = 0.D0 - fPsi = self%fPsi(xii) + fPsi = self%fPsi(Xi) r = (/DOT_PRODUCT(fPsi, self%x), & DOT_PRODUCT(fPsi, self%y), & DOT_PRODUCT(fPsi, self%z)/) @@ -228,17 +228,17 @@ MODULE moduleMesh3DCart END FUNCTION randPosEdgeTria !Shape functions for triangular surface - PURE FUNCTION fPsiEdgeTria(xii) RESULT(fPsi) + PURE FUNCTION fPsiEdgeTria(Xi) RESULT(fPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xii(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8), ALLOCATABLE:: fPsi(:) ALLOCATE(fPsi(1:3)) - fPsi(1) = 1.D0 - xii(1) - xii(2) - fPsi(2) = xii(1) - fPsi(3) = xii(2) + fPsi(1) = 1.D0 - Xi(1) - Xi(2) + fPsi(2) = Xi(1) + fPsi(3) = Xi(2) END FUNCTION fPsiEdgeTria @@ -254,6 +254,7 @@ MODULE moduleMesh3DCart INTEGER, INTENT(in):: p(:) TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 !Positions of each node + REAL(8):: Xi(1:3), fPsi(1:4) REAL(8):: volNodes(1:4) !Volume of each node self%n = n @@ -274,8 +275,9 @@ MODULE moduleMesh3DCart CALL self%calcVol() !Assign proportional volume to each node - !TODO: Review this to apply to all elements in the future - volNodes = self%fPsi((/0.25D0, 0.25D0, 0.25D0/))*self%volume + Xi = (/0.25D0, 0.25D0, 0.25D0/) + CALL self%fPsi(Xi, fPsi) + volNodes = fPsi*self%volume self%n1%v = self%n1%v + volNodes(1) self%n2%v = self%n2%v + volNodes(2) self%n3%v = self%n3%v + volNodes(3) @@ -295,19 +297,18 @@ MODULE moduleMesh3DCart CLASS(meshVol3DCartTetra), INTENT(in):: self REAL(8):: r(1:3) - REAL(8):: xii(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8):: Xi(1:3) + REAL(8):: fPsi(1:4) - xii(1) = random( 0.D0, 1.D0) - xii(2) = random( 0.D0, 1.D0 - xii(1)) - xii(3) = random( 0.D0, 1.D0 - xii(1) - xii(2)) + Xi(1) = random( 0.D0, 1.D0) + Xi(2) = random( 0.D0, 1.D0 - Xi(1)) + Xi(3) = random( 0.D0, 1.D0 - Xi(1) - Xi(2)) - ALLOCATE(fPsi(1:4)) - fPsi = self%fPsi(xii) + CALL self%fPsi(Xi, fPsi) - r(1) = DOT_PRODUCT(fPsi, self%x) - r(2) = DOT_PRODUCT(fPsi, self%y) - r(3) = DOT_PRODUCT(fPsi, self%z) + r = (/ DOT_PRODUCT(fPsi, self%x), & + DOT_PRODUCT(fPsi, self%y), & + DOT_PRODUCT(fPsi, self%z) /) END FUNCTION randPosVolTetra @@ -316,83 +317,81 @@ MODULE moduleMesh3DCart IMPLICIT NONE CLASS(meshVol3DCartTetra), INTENT(inout):: self - REAL(8):: xii(1:3) + REAL(8):: Xi(1:3) self%volume = 0.D0 - xii = (/0.25D0, 0.25D0, 0.25D0/) - self%volume = self%detJac(xii) + Xi = (/0.25D0, 0.25D0, 0.25D0/) + self%volume = self%detJac(Xi) END SUBROUTINE volumeTetra - !Computes element functions in point xii - PURE FUNCTION fPsiTetra(xi) RESULT(fPsi) + !Computes element functions in point Xi + PURE SUBROUTINE fPsiTetra(Xi, fPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8), INTENT(in):: Xi(1:3) + REAL(8), INTENT(out):: fPsi(:) - ALLOCATE(fPsi(1:4)) + fPsi(1) = 1.D0 - Xi(1) - Xi(2) - Xi(3) + fPsi(2) = Xi(1) + fPsi(3) = Xi(2) + fPsi(4) = Xi(3) - fPsi(1) = 1.D0 - xi(1) - xi(2) - xi(3) - fPsi(2) = xi(1) - fPsi(3) = xi(2) - fPsi(4) = xi(3) + END SUBROUTINE fPsiTetra - END FUNCTION fPsiTetra - - !Derivative element function at coordinates xii - PURE FUNCTION dPsiTetra(xii) RESULT(dPsi) + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiTetra(Xi) RESULT(dPsi) IMPLICIT NONE - REAL(8), INTENT(in):: xii(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8), ALLOCATABLE:: dPsi(:,:) ALLOCATE(dPsi(1:3,1:4)) - dPsi(1,:) = dPsiTetraXii1(xii(2), xii(3)) - dPsi(2,:) = dPsiTetraXii2(xii(1), xii(3)) - dPsi(3,:) = dPsiTetraXii3(xii(1), xii(2)) + dPsi(1,:) = dPsiTetraXi1(Xi(2), Xi(3)) + dPsi(2,:) = dPsiTetraXi2(Xi(1), Xi(3)) + dPsi(3,:) = dPsiTetraXi3(Xi(1), Xi(2)) END FUNCTION dPsiTetra - !Derivative element function respect to xii1 - PURE FUNCTION dPsiTetraXii1(xii2, xii3) RESULT(dPsiXii1) + !Derivative element function respect to Xi1 + PURE FUNCTION dPsiTetraXi1(Xi2, Xi3) RESULT(dPsiXi1) IMPLICIT NONE - REAL(8), INTENT(in):: xii2, xii3 - REAL(8):: dPsiXii1(1:4) + REAL(8), INTENT(in):: Xi2, Xi3 + REAL(8):: dPsiXi1(1:4) - dPsiXii1(1) = -1.D0 - dPsiXii1(2) = 1.D0 - dPsiXii1(3) = 0.D0 - dPsiXii1(4) = 0.D0 + dPsiXi1(1) = -1.D0 + dPsiXi1(2) = 1.D0 + dPsiXi1(3) = 0.D0 + dPsiXi1(4) = 0.D0 - END FUNCTION dPsiTetraXii1 + END FUNCTION dPsiTetraXi1 - !Derivative element function respect to xii2 - PURE FUNCTION dPsiTetraXii2(xii1, xii3) RESULT(dPsiXii2) + !Derivative element function respect to Xi2 + PURE FUNCTION dPsiTetraXi2(Xi1, Xi3) RESULT(dPsiXi2) IMPLICIT NONE - REAL(8), INTENT(in):: xii1, xii3 - REAL(8):: dPsiXii2(1:4) + REAL(8), INTENT(in):: Xi1, Xi3 + REAL(8):: dPsiXi2(1:4) - dPsiXii2(1) = -1.D0 - dPsiXii2(2) = 0.D0 - dPsiXii2(3) = 1.D0 - dPsiXii2(4) = 0.D0 + dPsiXi2(1) = -1.D0 + dPsiXi2(2) = 0.D0 + dPsiXi2(3) = 1.D0 + dPsiXi2(4) = 0.D0 - END FUNCTION dPsiTetraXii2 + END FUNCTION dPsiTetraXi2 - !Derivative element function respect to xii3 - PURE FUNCTION dPsiTetraXii3(xii1, xii2) RESULT(dPsiXii3) + !Derivative element function respect to Xi3 + PURE FUNCTION dPsiTetraXi3(Xi1, Xi2) RESULT(dPsiXi3) IMPLICIT NONE - REAL(8), INTENT(in):: xii1, xii2 - REAL(8):: dPsiXii3(1:4) + REAL(8), INTENT(in):: Xi1, Xi2 + REAL(8):: dPsiXi3(1:4) - dPsiXii3(1) = -1.D0 - dPsiXii3(2) = 0.D0 - dPsiXii3(3) = 0.D0 - dPsiXii3(4) = 1.D0 + dPsiXi3(1) = -1.D0 + dPsiXi3(2) = 0.D0 + dPsiXi3(3) = 0.D0 + dPsiXi3(4) = 1.D0 - END FUNCTION dPsiTetraXii3 + END FUNCTION dPsiTetraXi3 !Computes the derivatives in global coordinates PURE SUBROUTINE partialDerTetra(self, dPsi, dx, dy, dz) @@ -421,19 +420,19 @@ MODULE moduleMesh3DCart CLASS(meshVol3DCartTetra), INTENT(in):: self REAL(8), ALLOCATABLE:: localK(:,:) - REAL(8):: xii(1:3) + REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) REAL(8):: invJ(1:3,1:3), detJ ALLOCATE(localK(1:4,1:4)) localK = 0.D0 - xii = 0.D0 + Xi = 0.D0 !TODO: One point Gauss integral. Upgrade when possible - xii = (/ 0.25D0, 0.25D0, 0.25D0 /) - dPsi = self%dPsi(xii) - detJ = self%detJac(xii, dPsi) - invJ = self%invJac(xii, dPsi) - fPsi = self%fPsi(xii) + Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi, dPsi) + invJ = self%invJac(Xi, dPsi) + CALL self%fPsi(Xi, fPsi) localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ END FUNCTION elemKTetra @@ -445,40 +444,40 @@ MODULE moduleMesh3DCart REAL(8), INTENT(in):: source(1:) REAL(8), ALLOCATABLE:: localF(:) REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) - REAL(8):: xii(1:3) + REAL(8):: Xi(1:3) REAL(8):: detJ, f ALLOCATE(localF(1:4)) + localF = 0.D0 - xii = 0.D0 - !TODO: One point Gauss integral. Upgrade when possible - xii = (/ 0.25D0, 0.25D0, 0.25D0 /) - dPsi = self%dPsi(xii) - detJ = self%detJac(xii, dPsi) - fPsi = self%fPsi(xii) + Xi = 0.D0 + Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi, dPsi) + CALL self%fPsi(Xi, fPsi) f = DOT_PRODUCT(fPsi, source) localF = f*fPsi*1.D0*detJ END FUNCTION elemFTetra - PURE FUNCTION insideTetra(xi) RESULT(ins) + PURE FUNCTION insideTetra(Xi) RESULT(ins) IMPLICIT NONE - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) LOGICAL:: ins - ins = xi(1) >= 0.D0 .AND. & - xi(2) >= 0.D0 .AND. & - xi(3) >= 0.D0 .AND. & - 1.D0 - xi(1) - xi(2) - xi(3) >= 0.D0 + ins = Xi(1) >= 0.D0 .AND. & + Xi(2) >= 0.D0 .AND. & + Xi(3) >= 0.D0 .AND. & + 1.D0 - Xi(1) - Xi(2) - Xi(3) >= 0.D0 END FUNCTION insideTetra - PURE FUNCTION gatherEFTetra(self, xi) RESULT(EF) + PURE FUNCTION gatherEFTetra(self, Xi) RESULT(EF) IMPLICIT NONE CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8):: dPsi(1:3, 1:4) REAL(8):: dPsiR(1:3, 1:4) REAL(8):: invJ(1:3, 1:3), detJ @@ -490,9 +489,9 @@ MODULE moduleMesh3DCart self%n3%emData%phi, & self%n4%emData%phi /) - dPsi = self%dPsi(xi) - detJ = self%detJac(xi, dPsi) - invJ = self%invJac(xi, dPsi) + dPsi = self%dPsi(Xi) + detJ = self%detJac(Xi, dPsi) + invJ = self%invJac(Xi, dPsi) dPsiR = MATMUL(invJ, dPsi)/detJ EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi) EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi) @@ -500,11 +499,11 @@ MODULE moduleMesh3DCart END FUNCTION gatherEFTetra - PURE FUNCTION gatherMFTetra(self, xi) RESULT(MF) + PURE FUNCTION gatherMFTetra(self, Xi) RESULT(MF) IMPLICIT NONE CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) REAL(8):: fPsi(1:4) REAL(8):: MF_Nodes(1:4,1:3) REAL(8):: MF(1:3) @@ -522,7 +521,7 @@ MODULE moduleMesh3DCart self%n3%emData%B(3), & self%n4%emData%B(3) /) - fPsi = self%fPsi(xi) + CALL self%fPsi(Xi, fPsi) MF = MATMUL(fPsi, MF_Nodes) END FUNCTION gatherMFTetra @@ -538,37 +537,37 @@ MODULE moduleMesh3DCart END FUNCTION getNodesTetra - PURE FUNCTION phy2logTetra(self,r) RESULT(xi) + PURE FUNCTION phy2logTetra(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshVol3DCartTetra), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xi(1:3) + REAL(8):: Xi(1:3) REAL(8):: invJ(1:3, 1:3), detJ REAL(8):: deltaR(1:3) REAL(8):: dPsi(1:3, 1:4) - xi = 0.D0 + Xi = 0.D0 deltaR = (/r(1) - self%x(1), r(2) - self%y(1), r(3) - self%z(1) /) - dPsi = self%dPsi(xi) - invJ = self%invJac(xi, dPsi) - detJ = self%detJac(xi, dPsi) - xi = MATMUL(invJ, deltaR)/detJ + dPsi = self%dPsi(Xi) + invJ = self%invJac(Xi, dPsi) + detJ = self%detJac(Xi, dPsi) + Xi = MATMUL(invJ, deltaR)/detJ END FUNCTION phy2logTetra - SUBROUTINE nextElementTetra(self, xi, nextElement) + SUBROUTINE nextElementTetra(self, Xi, nextElement) IMPLICIT NONE CLASS(meshVol3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: xi(1:3) + REAL(8), INTENT(in):: Xi(1:3) CLASS(meshElement), POINTER, INTENT(out):: nextElement - REAL(8):: xiArray(1:4) + REAL(8):: XiArray(1:4) INTEGER:: nextInt !TODO: Review when connectivity - xiArray = (/ xi(3), 1.D0 - xi(1) - xi(2) - xi(3), xi(2), xi(1) /) - nextInt = MINLOC(xiArray, 1) + XiArray = (/ Xi(3), 1.D0 - Xi(1) - Xi(2) - Xi(3), Xi(2), Xi(1) /) + nextInt = MINLOC(XiArray, 1) NULLIFY(nextElement) SELECT CASE(nextInt) CASE (1) @@ -585,11 +584,11 @@ MODULE moduleMesh3DCart !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 3D !Computes element Jacobian determinant - PURE FUNCTION detJ3DCart(self, xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ3DCart(self, Xi, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshVol3DCart), 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:, 1:) REAL(8):: dJ REAL(8), ALLOCATABLE:: dPsi(:,:) @@ -599,7 +598,7 @@ MODULE moduleMesh3DCart dPsi = dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi) END IF @@ -610,11 +609,11 @@ MODULE moduleMesh3DCart END FUNCTION detJ3DCart - PURE FUNCTION invJ3DCart(self,xi,dPsi_in) RESULT(invJ) + PURE FUNCTION invJ3DCart(self,Xi,dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshVol3DCart), 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:,1:) REAL(8), ALLOCATABLE:: dPsi(:,:) REAL(8), DIMENSION(1:3):: dx, dy, dz @@ -624,7 +623,7 @@ MODULE moduleMesh3DCart dPsi=dPsi_in ELSE - dPsi = self%dPsi(xi) + dPsi = self%dPsi(Xi) END IF diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 05c4b85..11e8bc7 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -211,11 +211,11 @@ MODULE moduleMesh END FUNCTION getNodesVol_interface - PURE FUNCTION fPsi_interface(xi) RESULT(fPsi) + PURE SUBROUTINE fPsi_interface(xi, fPsi) REAL(8), INTENT(in):: xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8), INTENT(out):: fPsi(:) - END FUNCTION fPsi_interface + END SUBROUTINE fPsi_interface PURE FUNCTION elemK_interface(self) RESULT(localK) IMPORT:: meshVol @@ -496,11 +496,14 @@ MODULE moduleMesh INTEGER:: i, nNodes CLASS(meshNode), POINTER:: node - fPsi = self%fPsi(part%xi) - tensorS = outerProduct(part%v, part%v) - sp = part%species%n volNodes = self%getNodes() nNodes = SIZE(volNodes) + ALLOCATE(fPsi(1:nNodes)) + CALL self%fPsi(part%xi, fPsi) + + tensorS = outerProduct(part%v, part%v) + + sp = part%species%n DO i = 1, nNodes node => mesh%nodes(volNodes(i))%obj