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