diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 1a439cc..0f409fb 100644 --- a/src/modules/init/moduleInput.f90 +++ b/src/modules/init/moduleInput.f90 @@ -359,8 +359,8 @@ MODULE moduleInput DO e = 1, mesh%numCells !Scale variables !Density at centroid of cell - nodes = mesh%cells(e)%obj%getNodes() nNodes = mesh%cells(e)%obj%nNodes + nodes = mesh%cells(e)%obj%getNodes(nNodes) ALLOCATE(fPsi(1:nNodes)) fPsi = mesh%cells(e)%obj%fPsi((/0.D0, 0.D0, 0.D0/), nNodes) ALLOCATE(source(1:nNodes)) diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index 65b13a3..5f89f20 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -89,11 +89,12 @@ MODULE moduleMesh0D END SUBROUTINE initCell0D - PURE FUNCTION getNodes0D(self) RESULT(n) + PURE FUNCTION getNodes0D(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell0D), INTENT(in):: self - INTEGER:: n(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) n = self%n1%n @@ -133,46 +134,50 @@ MODULE moduleMesh0D END FUNCTION dPsi0D - PURE FUNCTION detJ0D(self, Xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ0D(self, Xi, nNodes, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshCell0D), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) REAL(8):: dJ dJ = 0.D0 END FUNCTION detJ0D - PURE FUNCTION invJ0D(self, Xi, dPsi_in) RESULT(invJ) + PURE FUNCTION invJ0D(self, Xi, nNodes, dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshCell0D), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) REAL(8):: invJ(1:3,1:3) invJ = 0.D0 END FUNCTION invJ0D - PURE FUNCTION elemK0D(self) RESULT(localK) + PURE FUNCTION elemK0D(self, nNodes) RESULT(localK) IMPLICIT NONE CLASS(meshCell0D), INTENT(in):: self - REAL(8):: localK(1:self%nNodes,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) localK = 0.D0 END FUNCTION elemK0D - PURE FUNCTION elemF0D(self, source) RESULT(localF) + PURE FUNCTION elemF0D(self, nNodes, source) RESULT(localF) IMPLICIT NONE CLASS(meshCell0D), INTENT(in):: self - REAL(8), INTENT(in):: source(1:self%nNodes) - REAL(8):: localF(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) localF = 0.D0 @@ -187,7 +192,7 @@ MODULE moduleMesh0D phi = (/ self%n1%emData%phi /) - array = -self%gatherDF(Xi, phi) + array = -self%gatherDF(Xi, 1, phi) END FUNCTION gatherEF0D @@ -204,7 +209,7 @@ MODULE moduleMesh0D B(:,3) = (/ self%n1%emData%B(3) /) - array = self%gatherF(Xi, 3, B) + array = self%gatherF(Xi, 1, B) END FUNCTION gatherMF0D diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 982c703..56d1b65 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -41,10 +41,11 @@ MODULE moduleMesh1DCart END TYPE meshCell1DCart ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, dPsi, dx) + PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx) IMPORT meshCell1DCart CLASS(meshCell1DCart), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1):: dx END SUBROUTINE partialDer_interface @@ -129,6 +130,7 @@ MODULE moduleMesh1DCart INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj !Get element coordinates r1 = self%n1%getCoordinates() @@ -152,13 +154,13 @@ MODULE moduleMesh1DCart END SUBROUTINE initEdge1DCart !Get nodes from edge - PURE FUNCTION getNodes1DCart(self) RESULT(n) + PURE FUNCTION getNodes1DCart(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge1DCart), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1)) n = (/ self%n1%n /) END FUNCTION getNodes1DCart @@ -250,7 +252,7 @@ MODULE moduleMesh1DCart !1 point Gauss integral Xi = 0.D0 fPsi = self%fPsi(Xi, 2) - detJ = self%detJac(Xi) + detJ = self%detJac(Xi, 2) l = 2.D0*detJ self%volume = l self%arNodes = fPsi*l @@ -290,11 +292,12 @@ MODULE moduleMesh1DCart END FUNCTION dPsiSegm !Computes partial derivatives of coordinates - PURE SUBROUTINE partialDerSegm(self, dPsi, dx) + PURE SUBROUTINE partialDerSegm(self, nNodes, dPsi, dx) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1):: dx dx(1) = DOT_PRODUCT(dPsi(1,:), self%x) @@ -302,11 +305,12 @@ MODULE moduleMesh1DCart END SUBROUTINE partialDerSegm !Computes local stiffness matrix - PURE FUNCTION elemKSegm(self) RESULT(localK) + PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self - REAL(8):: localK(1:self%nNodes,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) REAL(8):: Xi(1:3) REAL(8):: dPsi(1:3, 1:2) REAL(8):: invJ(1:3,1:3), detJ @@ -317,8 +321,8 @@ MODULE moduleMesh1DCart DO l = 1, 3 Xi(1) = corSeg(l) dPsi = self%dPsi(Xi, 2) - detJ = self%detJac(Xi, dPsi) - invJ = self%invJac(Xi, dPsi) + detJ = self%detJac(Xi, 2, dPsi) + invJ = self%invJac(Xi, 2, dPsi) localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & wSeg(l)/detJ @@ -327,12 +331,13 @@ MODULE moduleMesh1DCart END FUNCTION elemKSegm - PURE FUNCTION elemFSegm(self, source) RESULT(localF) + PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self - REAL(8), INTENT(in):: source(1:self%nNodes) - REAL(8):: localF(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: fPsi(1:2) REAL(8):: detJ, f REAL(8):: Xi(1:3) @@ -343,7 +348,7 @@ MODULE moduleMesh1DCart DO l = 1, 3 Xi(1) = corSeg(l) - detJ = self%detJac(Xi) + detJ = self%detJac(Xi, 2) fPsi = self%fPsi(Xi, 2) f = DOT_PRODUCT(fPsi, source) localF = localF + f*fPsi*wSeg(l)*detJ @@ -362,7 +367,7 @@ MODULE moduleMesh1DCart phi = (/ self%n1%emData%phi, & self%n2%emData%phi /) - array = -self%gatherDF(Xi, phi) + array = -self%gatherDF(Xi, 2, phi) END FUNCTION gatherEFSegm @@ -382,7 +387,7 @@ MODULE moduleMesh1DCart B(:,3) = (/ self%n1%emData%B(3), & self%n2%emData%B(3) /) - array = self%gatherF(Xi, 3, B) + array = self%gatherF(Xi, 2, B) END FUNCTION gatherMFSegm @@ -398,11 +403,12 @@ MODULE moduleMesh1DCart END FUNCTION insideSegm !Get nodes from 1D volume - PURE FUNCTION getNodesSegm(self) RESULT(n) + PURE FUNCTION getNodesSegm(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self - INTEGER:: n(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) n = (/ self%n1%n, self%n2%n /) @@ -442,13 +448,14 @@ MODULE moduleMesh1DCart !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS !Calculates a random position in 1D volume !Computes the element Jacobian determinant - PURE FUNCTION detJ1DCart(self, Xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ1DCart(self, Xi, nNodes, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshCell1DCart), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) - REAL(8):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) + REAL(8):: dPsi(1:3,1:nNodes) REAL(8):: dJ REAL(8):: dx(1) @@ -460,20 +467,21 @@ MODULE moduleMesh1DCart END IF - CALL self%partialDer(dPsi, dx) + CALL self%partialDer(2, dPsi, dx) dJ = dx(1) END FUNCTION detJ1DCart !Computes the invers Jacobian - PURE FUNCTION invJ1DCart(self, Xi, dPsi_in) RESULT(invJ) + PURE FUNCTION invJ1DCart(self, Xi, nNodes, dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshCell1DCart), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) REAL(8):: invJ(1:3,1:3) - REAL(8):: dPsi(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:nNodes) REAL(8):: dx(1) IF (PRESENT(dPsi_in)) THEN @@ -486,7 +494,7 @@ MODULE moduleMesh1DCart invJ = 0.D0 - CALL self%partialDer(dPsi, dx) + CALL self%partialDer(2, dPsi, dx) invJ(1,1) = 1.D0/dx(1) diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 1ce7836..8b441be 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -41,11 +41,12 @@ MODULE moduleMesh1DRad END TYPE meshCell1DRad ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, dPsi, dx) + PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx) IMPORT meshCell1DRad CLASS(meshCell1DRad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1):: dx END SUBROUTINE partialDer_interface @@ -130,6 +131,7 @@ MODULE moduleMesh1DRad INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj !Get element coordinates r1 = self%n1%getCoordinates() @@ -154,13 +156,13 @@ MODULE moduleMesh1DRad END SUBROUTINE initEdge1DRad !Get nodes from edge - PURE FUNCTION getNodes1DRad(self) RESULT(n) + PURE FUNCTION getNodes1DRad(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge1DRad), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1)) n = (/ self%n1%n /) END FUNCTION getNodes1DRad @@ -253,17 +255,17 @@ MODULE moduleMesh1DRad !1 point Gauss integral Xi = 0.D0 fPsi = self%fPsi(Xi, 2) - detJ = self%detJac(Xi) + detJ = self%detJac(Xi, 2) !Computes total volume of the cell r = DOT_PRODUCT(fPsi, self%r) l = 2.D0*detJ self%volume = r*l !Computes volume per node Xi = (/-5.D-1, 0.D0, 0.D0/) - r = self%gatherF(Xi, self%r) + r = self%gatherF(Xi, 2, self%r) self%arNodes(1) = fPsi(1)*r*l Xi = (/ 5.D-1, 0.D0, 0.D0/) - r = self%gatherF(Xi, self%r) + r = self%gatherF(Xi, 2, self%r) self%arNodes(2) = fPsi(2)*r*l END SUBROUTINE areaRad @@ -301,11 +303,12 @@ MODULE moduleMesh1DRad END FUNCTION dPsiRad !Computes partial derivatives of coordinates - PURE SUBROUTINE partialDerRad(self, dPsi, dx) + PURE SUBROUTINE partialDerRad(self, nNodes, dPsi, dx) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1):: dx dx(1) = DOT_PRODUCT(dPsi(1,:), self%r) @@ -313,12 +316,13 @@ MODULE moduleMesh1DRad END SUBROUTINE partialDerRad !Computes local stiffness matrix - PURE FUNCTION elemKRad(self) RESULT(localK) + PURE FUNCTION elemKRad(self, nNodes) RESULT(localK) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self - REAL(8):: localK(1:self%nNodes,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) REAL(8):: Xi(1:3) REAL(8):: dPsi(1:3, 1:2) REAL(8):: invJ(1:3,1:3), detJ @@ -330,8 +334,8 @@ MODULE moduleMesh1DRad DO l = 1, 3 Xi(1) = corSeg(l) dPsi = self%dPsi(Xi, 2) - detJ = self%detJac(Xi, dPsi) - invJ = self%invJac(Xi, dPsi) + detJ = self%detJac(Xi, 2, dPsi) + invJ = self%invJac(Xi, 2, dPsi) fPsi = self%fPsi(Xi, 2) r = DOT_PRODUCT(fPsi, self%r) localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & @@ -344,13 +348,14 @@ MODULE moduleMesh1DRad END FUNCTION elemKRad - PURE FUNCTION elemFRad(self, source) RESULT(localF) + PURE FUNCTION elemFRad(self, nNodes, source) RESULT(localF) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self - REAL(8), INTENT(in):: source(1:self%nNodes) - REAL(8):: localF(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: fPsi(1:2) REAL(8):: detJ, f, r REAL(8):: Xi(1:3) @@ -361,7 +366,7 @@ MODULE moduleMesh1DRad DO l = 1, 3 Xi(1) = corSeg(l) - detJ = self%detJac(Xi) + detJ = self%detJac(Xi, 2) fPsi = self%fPsi(Xi, 2) r = DOT_PRODUCT(fPsi, self%r) f = DOT_PRODUCT(fPsi, source) @@ -381,7 +386,7 @@ MODULE moduleMesh1DRad phi = (/ self%n1%emData%phi, & self%n2%emData%phi /) - array = -self%gatherDF(Xi, phi) + array = -self%gatherDF(Xi, 2, phi) END FUNCTION gatherEFRad @@ -401,7 +406,7 @@ MODULE moduleMesh1DRad B(:,3) = (/ self%n1%emData%B(3), & self%n2%emData%B(3) /) - array = self%gatherF(Xi, 3, B) + array = self%gatherF(Xi, 2, B) END FUNCTION gatherMFRad @@ -417,11 +422,12 @@ MODULE moduleMesh1DRad END FUNCTION insideRad !Get nodes from 1D volume - PURE FUNCTION getNodesRad(self) RESULT(n) + PURE FUNCTION getNodesRad(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self - INTEGER:: n(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) n = (/ self%n1%n, self%n2%n /) @@ -460,13 +466,14 @@ MODULE moduleMesh1DRad !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS !Computes the element Jacobian determinant - PURE FUNCTION detJ1DRad(self, Xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ1DRad(self, Xi, nNodes, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshCell1DRad), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) - REAL(8):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) + REAL(8):: dPsi(1:3,1:nNodes) REAL(8):: dJ REAL(8):: dx(1) @@ -478,19 +485,20 @@ MODULE moduleMesh1DRad END IF - CALL self%partialDer(dPsi, dx) + CALL self%partialDer(nNodes, dPsi, dx) dJ = dx(1) END FUNCTION detJ1DRad !Computes the invers Jacobian - PURE FUNCTION invJ1DRad(self, Xi, dPsi_in) RESULT(invJ) + PURE FUNCTION invJ1DRad(self, Xi, nNodes, dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshCell1DRad), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) - REAL(8):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) + REAL(8):: dPsi(1:3,1:nNodes) REAL(8):: dx(1) REAL(8):: invJ(1:3,1:3) @@ -504,7 +512,7 @@ MODULE moduleMesh1DRad invJ = 0.D0 - CALL self%partialDer(dPsi, dx) + CALL self%partialDer(nNodes, dPsi, dx) invJ(1,1) = 1.D0/dx(1) diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 649b8fb..2ccde7d 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -46,10 +46,11 @@ MODULE moduleMesh2DCart END TYPE meshCell2DCart ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, dPsi, dx, dy) + PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx, dy) IMPORT meshCell2DCart CLASS(meshCell2DCart), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy END SUBROUTINE partialDer_interface @@ -166,6 +167,7 @@ MODULE moduleMesh2DCart INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj self%n2 => mesh%nodes(p(2))%obj !Get element coordinates @@ -194,13 +196,13 @@ MODULE moduleMesh2DCart END SUBROUTINE initEdge2DCart !Get nodes from edge - PURE FUNCTION getNodes2DCart(self) RESULT(n) + PURE FUNCTION getNodes2DCart(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge2DCart), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1:2)) n = (/self%n1%n, self%n2%n /) END FUNCTION getNodes2DCart @@ -255,8 +257,13 @@ MODULE moduleMesh2DCart TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:) REAL(8), DIMENSION(1:3):: r1, r2, r3, r4 + !Assign node index self%n = n + + !Assign number of nodes of cell self%nNodes = SIZE(p) + + !Assign nodes to element self%n1 => nodes(p(1))%obj self%n2 => nodes(p(2))%obj self%n3 => nodes(p(3))%obj @@ -296,7 +303,7 @@ MODULE moduleMesh2DCart self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral Xi = 0.D0 - detJ = self%detJac(Xi)*4.D0 !4 + detJ = self%detJac(Xi, 4)*4.D0 !4 fPsi = self%fPsi(Xi, 4) self%volume = detJ self%arNodes = fPsi*detJ @@ -347,11 +354,12 @@ MODULE moduleMesh2DCart END FUNCTION dPsiQuad !Partial derivative in global coordinates - PURE SUBROUTINE partialDerQuad(self, dPsi, dx, dy) + PURE SUBROUTINE partialDerQuad(self, nNodes, dPsi, dx, dy) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy dx = (/ DOT_PRODUCT(dPsi(1,1:4),self%x(1:4)), & @@ -384,11 +392,12 @@ MODULE moduleMesh2DCart END FUNCTION randPosCellQuad !Computes element local stiffness matrix - PURE FUNCTION elemKQuad(self) RESULT(localK) + PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self - REAL(8):: localK(1:self%nNodes,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:3,1:4) REAL(8):: invJ(1:3,1:3), detJ @@ -403,8 +412,8 @@ MODULE moduleMesh2DCart Xi(1) = corQuad(m) fPsi = self%fPsi(Xi, 4) dPsi = self%dPsi(Xi, 4) - detJ = self%detJac(Xi,dPsi) - invJ = self%invJac(Xi,dPsi) + detJ = self%detJac(Xi, 4, dPsi) + invJ = self%invJac(Xi, 4, dPsi) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & MATMUL(invJ,dPsi))* & wQuad(l)*wQuad(m)/detJ @@ -415,12 +424,13 @@ MODULE moduleMesh2DCart END FUNCTION elemKQuad !Computes the local source vector for a force f - PURE FUNCTION elemFQuad(self, source) RESULT(localF) + PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: source(1:self%nNodes) - REAL(8):: localF(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) REAL(8):: detJ, f @@ -432,7 +442,7 @@ MODULE moduleMesh2DCart Xi(1) = corQuad(l) DO m = 1, 3 Xi(2) = corQuad(m) - detJ = self%detJac(Xi) + detJ = self%detJac(Xi, 4) fPsi = self%fPsi(Xi, 4) f = DOT_PRODUCT(fPsi,source) localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ @@ -454,7 +464,7 @@ MODULE moduleMesh2DCart self%n3%emData%phi, & self%n4%emData%phi /) - array = -self%gatherDF(Xi, phi) + array = -self%gatherDF(Xi, 4, phi) END FUNCTION gatherEFQuad @@ -480,7 +490,7 @@ MODULE moduleMesh2DCart self%n3%emData%B(3), & self%n4%emData%B(3) /) - array = self%gatherF(Xi, 3, B) + array = self%gatherF(Xi, 4, B) END FUNCTION gatherMFQuad @@ -497,11 +507,12 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad !Gets nodes from quadrilateral element - PURE FUNCTION getNodesQuad(self) RESULT(n) + PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self - INTEGER:: n(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) @@ -524,7 +535,7 @@ MODULE moduleMesh2DCart DO WHILE(conv > 1.D-2) dPsi = self%dPsi(XiO, 4) - invJ = self%invJac(XiO, dPsi) + invJ = self%invJac(XiO, 4, dPsi) fPsi = self%fPsi(XiO, 4) f = (/ DOT_PRODUCT(fPsi,self%x), & DOT_PRODUCT(fPsi,self%y), & @@ -641,7 +652,7 @@ MODULE moduleMesh2DCart 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 + detJ = self%detJac(Xi, 4)/2.D0 fPsi = self%fPsi(Xi, 4) self%volume = detJ self%arNodes = fPsi*detJ @@ -679,11 +690,12 @@ MODULE moduleMesh2DCart END FUNCTION dPsiTria - PURE SUBROUTINE partialDerTria(self, dPsi, dx, dy) + PURE SUBROUTINE partialDerTria(self, nNodes, dPsi, dx, dy) IMPLICIT NONE CLASS(meshCell2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy dx = (/ DOT_PRODUCT(dPsi(1,:),self%x), & @@ -694,11 +706,12 @@ MODULE moduleMesh2DCart END SUBROUTINE partialDerTria !Computes element local stiffness matrix - PURE FUNCTION elemKTria(self) RESULT(localK) + PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) IMPLICIT NONE CLASS(meshCell2DCartTria), INTENT(in):: self - REAL(8):: localK(1:self%nNodes,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:3), dPsi(1:3,1:3) REAL(8):: invJ(1:3,1:3), detJ @@ -710,10 +723,10 @@ MODULE moduleMesh2DCart DO l=1, 4 Xi(1) = Xi1Tria(l) Xi(2) = Xi2Tria(l) - dPsi = self%dPsi(Xi, 4) - detJ = self%detJac(Xi,dPsi) - invJ = self%invJac(Xi,dPsi) - fPsi = self%fPsi(Xi, 4) + dPsi = self%dPsi(Xi, 3) + detJ = self%detJac(Xi, 3, dPsi) + invJ = self%invJac(Xi, 3, dPsi) + fPsi = self%fPsi(Xi, 3) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ END DO @@ -721,12 +734,13 @@ MODULE moduleMesh2DCart END FUNCTION elemKTria !Computes element local source vector - PURE FUNCTION elemFTria(self, source) RESULT(localF) + PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF) IMPLICIT NONE CLASS(meshCell2DCartTria), INTENT(in):: self - REAL(8), INTENT(in):: source(1:self%nNodes) - REAL(8):: localF(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: fPsi(1:3) REAL(8):: Xi(1:3) REAL(8):: detJ, f @@ -738,8 +752,8 @@ MODULE moduleMesh2DCart DO l=1, 4 Xi(1) = Xi1Tria(l) Xi(2) = Xi2Tria(l) - detJ = self%detJac(Xi) - fPsi = self%fPsi(Xi, 4) + detJ = self%detJac(Xi, 3) + fPsi = self%fPsi(Xi, 3) f = DOT_PRODUCT(fPsi,source) localF = localF + f*fPsi*wTria(l)*detJ @@ -758,7 +772,7 @@ MODULE moduleMesh2DCart self%n2%emData%phi, & self%n3%emData%phi /) - array = -self%gatherDF(Xi, phi) + array = -self%gatherDF(Xi, 3, phi) END FUNCTION gatherEFTria @@ -798,11 +812,12 @@ MODULE moduleMesh2DCart END FUNCTION insideTria !Gets node indexes from triangular element - PURE FUNCTION getNodesTria(self) RESULT(n) + PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell2DCartTria), INTENT(in):: self - INTEGER:: n(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) n = (/self%n1%n, self%n2%n, self%n3%n /) @@ -822,9 +837,9 @@ MODULE moduleMesh2DCart !Direct method to convert coordinates Xi = 0.D0 deltaR = (/ r(1) - self%x(1), r(2) - self%y(1), 0.D0 /) - dPsi = self%dPsi(Xi, 4) - invJ = self%invJac(Xi, dPsi) - detJ = self%detJac(Xi, dPsi) + dPsi = self%dPsi(Xi, 3) + invJ = self%invJac(Xi, 3, dPsi) + detJ = self%detJac(Xi, 3, dPsi) Xi = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria @@ -854,14 +869,15 @@ 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, nNodes, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshCell2DCart), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) REAL(8):: dJ - REAL(8):: dPsi(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:nNodes) REAL(8):: dx(1:2), dy(1:2) IF(PRESENT(dPsi_in)) THEN @@ -872,21 +888,22 @@ MODULE moduleMesh2DCart END IF - CALL self%partialDer(dPsi, dx, dy) + CALL self%partialDer(nNodes, dPsi, dx, dy) dJ = dx(1)*dy(2)-dx(2)*dy(1) END FUNCTION detJ2DCart !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJ2DCart(self,Xi,dPsi_in) RESULT(invJ) + PURE FUNCTION invJ2DCart(self, Xi, nNodes, dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshCell2DCart), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) REAL(8):: invJ(1:3,1:3) - REAL(8):: dPsi(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:nNodes) REAL(8):: dx(1:2), dy(1:2) IF(PRESENT(dPsi_in)) THEN @@ -899,7 +916,7 @@ MODULE moduleMesh2DCart invJ = 0.D0 - CALL self%partialDer(dPsi, dx, dy) + CALL self%partialDer(nNodes, dPsi, dx, dy) invJ(1,1:2) = (/ dy(2), -dx(2) /) invJ(2,1:2) = (/ -dy(1), dx(1) /) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index adf5d4f..d33fbbf 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -46,10 +46,11 @@ MODULE moduleMesh2DCyl END TYPE meshCell2DCyl ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, dPsi, dz, dr) + PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dz, dr) IMPORT meshCell2DCyl CLASS(meshCell2DCyl), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr END SUBROUTINE partialDer_interface @@ -166,6 +167,7 @@ MODULE moduleMesh2DCyl INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj self%n2 => mesh%nodes(p(2))%obj !Get element coordinates @@ -195,13 +197,13 @@ MODULE moduleMesh2DCyl END SUBROUTINE initEdge2DCyl !Get nodes from edge - PURE FUNCTION getNodes2DCyl(self) RESULT(n) + PURE FUNCTION getNodes2DCyl(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge2DCyl), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1:2)) n = (/self%n1%n, self%n2%n /) END FUNCTION getNodes2DCyl @@ -306,23 +308,23 @@ MODULE moduleMesh2DCyl self%arNodes = 0.D0 !2D 1 point Gauss Quad Integral Xi = 0.D0 - detJ = self%detJac(Xi)*PI8 !4*2*pi + detJ = self%detJac(Xi, 4)*PI8 !4*2*pi fPsi = self%fPsi(Xi, 4) !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 = self%gatherF(Xi, self%r) + r = self%gatherF(Xi, 4, self%r) self%arNodes(1) = fPsi(1)*r*detJ Xi = (/ 5.D-1, -5.D-1, 0.D0/) - r = self%gatherF(Xi, self%r) + r = self%gatherF(Xi, 4, self%r) self%arNodes(2) = fPsi(2)*r*detJ Xi = (/ 5.D-1, 5.D-1, 0.D0/) - r = self%gatherF(Xi, self%r) + r = self%gatherF(Xi, 4, self%r) self%arNodes(3) = fPsi(3)*r*detJ Xi = (/-5.D-1, 5.D-1, 0.D0/) - r = self%gatherF(Xi, self%r) + r = self%gatherF(Xi, 4, self%r) self%arNodes(4) = fPsi(4)*r*detJ END SUBROUTINE areaQuad @@ -371,11 +373,12 @@ MODULE moduleMesh2DCyl END FUNCTION dPsiQuad !Partial derivative in global coordinates - PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr) + PURE SUBROUTINE partialDerQuad(self, nNodes, dPsi, dz, dr) IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr dz = (/ DOT_PRODUCT(dPsi(1,1:4),self%z(1:4)), & @@ -408,12 +411,13 @@ MODULE moduleMesh2DCyl END FUNCTION randPosCellQuad !Computes element local stiffness matrix - PURE FUNCTION elemKQuad(self) RESULT(localK) + PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(in):: self - REAL(8):: localK(1:self%nNodes,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:3,1:4) REAL(8):: r @@ -429,8 +433,8 @@ MODULE moduleMesh2DCyl Xi(1) = corQuad(m) fPsi = self%fPsi(Xi, 4) dPsi = self%dPsi(Xi, 4) - detJ = self%detJac(Xi,dPsi) - invJ = self%invJac(Xi,dPsi) + detJ = self%detJac(Xi, 4, dPsi) + invJ = self%invJac(Xi, 4, dPsi) r = DOT_PRODUCT(fPsi,self%r) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & MATMUL(invJ,dPsi))* & @@ -443,13 +447,14 @@ MODULE moduleMesh2DCyl END FUNCTION elemKQuad !Computes the local source vector for a force f - PURE FUNCTION elemFQuad(self, source) RESULT(localF) + PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(in):: self - REAL(8), INTENT(in):: source(1:self%nNodes) - REAL(8):: localF(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) REAL(8):: r @@ -462,7 +467,7 @@ MODULE moduleMesh2DCyl Xi(1) = corQuad(l) DO m = 1, 3 Xi(2) = corQuad(m) - detJ = self%detJac(Xi) + detJ = self%detJac(Xi, 4) fPsi = self%fPsi(Xi, 4) r = DOT_PRODUCT(fPsi,self%r) f = DOT_PRODUCT(fPsi,source) @@ -486,7 +491,7 @@ MODULE moduleMesh2DCyl self%n3%emData%phi, & self%n4%emData%phi /) - array = -self%gatherDF(Xi, phi) + array = -self%gatherDF(Xi, 4, phi) END FUNCTION gatherEFQuad @@ -512,7 +517,7 @@ MODULE moduleMesh2DCyl self%n3%emData%B(3), & self%n4%emData%B(3) /) - array = self%gatherF(Xi, 3, B) + array = self%gatherF(Xi, 4, B) END FUNCTION gatherMFQuad @@ -529,11 +534,12 @@ MODULE moduleMesh2DCyl END FUNCTION insideQuad !Gets nodes from quadrilateral element - PURE FUNCTION getNodesQuad(self) RESULT(n) + PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(in):: self - INTEGER:: n(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) @@ -556,8 +562,8 @@ MODULE moduleMesh2DCyl DO WHILE(conv > 1.D-2) dPsi = self%dPsi(XiO, 4) - invJ = self%invJac(XiO, dPsi) - detJ = self%detJac(XiO, dPsi) + invJ = self%invJac(XiO, 4, dPsi) + detJ = self%detJac(XiO, 4, dPsi) fPsi = self%fPsi(XiO, 4) f = (/ DOT_PRODUCT(fPsi,self%z), & DOT_PRODUCT(fPsi,self%r), & @@ -676,7 +682,7 @@ MODULE moduleMesh2DCyl 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)*PI !2PI*1/2 + detJ = self%detJac(Xi, 3)*PI !2PI*1/2 fPsi = self%fPsi(Xi, 4) !Computes total volume of the cell r = DOT_PRODUCT(fPsi,self%r) @@ -717,11 +723,12 @@ MODULE moduleMesh2DCyl END FUNCTION dPsiTria - PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr) + PURE SUBROUTINE partialDerTria(self, nNodes, dPsi, dz, dr) IMPLICIT NONE CLASS(meshCell2DCylTria), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr dz = (/ DOT_PRODUCT(dPsi(1,:),self%z), & @@ -732,12 +739,13 @@ MODULE moduleMesh2DCyl END SUBROUTINE partialDerTria !Computes element local stiffness matrix - PURE FUNCTION elemKTria(self) RESULT(localK) + PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshCell2DCylTria), INTENT(in):: self - REAL(8):: localK(1:self%nNodes,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) REAL(8):: Xi(1:3) REAL(8):: r REAL(8):: fPsi(1:3), dPsi(1:3,1:3) @@ -751,8 +759,8 @@ MODULE moduleMesh2DCyl Xi(1) = Xi1Tria(l) Xi(2) = Xi2Tria(l) dPsi = self%dPsi(Xi, 4) - detJ = self%detJac(Xi,dPsi) - invJ = self%invJac(Xi,dPsi) + detJ = self%detJac(Xi, 3, dPsi) + invJ = self%invJac(Xi, 3, dPsi) fPsi = self%fPsi(Xi, 4) r = DOT_PRODUCT(fPsi,self%r) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ @@ -763,13 +771,14 @@ MODULE moduleMesh2DCyl END FUNCTION elemKTria !Computes element local source vector - PURE FUNCTION elemFTria(self, source) RESULT(localF) + PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF) USE moduleConstParam, ONLY: PI2 IMPLICIT NONE CLASS(meshCell2DCylTria), INTENT(in):: self - REAL(8), INTENT(in):: source(1:self%nNodes) - REAL(8):: localF(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: fPsi(1:3) REAL(8):: Xi(1:3) REAL(8):: r @@ -782,8 +791,8 @@ MODULE moduleMesh2DCyl DO l=1, 4 Xi(1) = Xi1Tria(l) Xi(2) = Xi2Tria(l) - detJ = self%detJac(Xi) - fPsi = self%fPsi(Xi, 4) + detJ = self%detJac(Xi, 3) + fPsi = self%fPsi(Xi, 3) r = DOT_PRODUCT(fPsi,self%r) f = DOT_PRODUCT(fPsi,source) localF = localF + r*f*fPsi*wTria(l)*detJ @@ -804,7 +813,7 @@ MODULE moduleMesh2DCyl self%n2%emData%phi, & self%n3%emData%phi /) - array = -self%gatherDF(Xi, phi) + array = -self%gatherDF(Xi, 4, phi) END FUNCTION gatherEFTria @@ -844,11 +853,12 @@ MODULE moduleMesh2DCyl END FUNCTION insideTria !Gets node indexes from triangular element - PURE FUNCTION getNodesTria(self) RESULT(n) + PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell2DCylTria), INTENT(in):: self - INTEGER:: n(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) n = (/self%n1%n, self%n2%n, self%n3%n /) @@ -868,9 +878,9 @@ MODULE moduleMesh2DCyl !Direct method to convert coordinates Xi = 0.D0 deltaR = (/ r(1) - self%z(1), r(2) - self%r(1), 0.D0 /) - dPsi = self%dPsi(Xi, 4) - invJ = self%invJac(Xi, dPsi) - detJ = self%detJac(Xi, dPsi) + dPsi = self%dPsi(Xi, 3) + invJ = self%invJac(Xi, 3, dPsi) + detJ = self%detJac(Xi, 3, dPsi) Xi = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria @@ -900,39 +910,41 @@ MODULE moduleMesh2DCyl !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS !Computes element Jacobian determinant - PURE FUNCTION detJ2DCyl(self, Xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ2DCyl(self, Xi, nNodes, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshCell2DCyl), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) REAL(8):: dJ - REAL(8):: dPsi(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:nNodes) REAL(8):: dz(1:2), dr(1:2) IF(PRESENT(dPsi_in)) THEN dPsi = dPsi_in ELSE - dPsi = self%dPsi(Xi, 4) + dPsi = self%dPsi(Xi, nNodes) END IF - CALL self%partialDer(dPsi, dz, dr) + CALL self%partialDer(nNodes, dPsi, dz, dr) dJ = dz(1)*dr(2)-dz(2)*dr(1) END FUNCTION detJ2DCyl !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJ2DCyl(self,Xi,dPsi_in) RESULT(invJ) + PURE FUNCTION invJ2DCyl(self, Xi, nNodes, dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshCell2DCyl), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) REAL(8):: invJ(1:3,1:3) - REAL(8):: dPsi(1:3,1:self%nNodes) + REAL(8):: dPsi(1:3,1:nNodes) REAL(8):: dz(1:2), dr(1:2) IF(PRESENT(dPsi_in)) THEN @@ -945,7 +957,7 @@ MODULE moduleMesh2DCyl invJ = 0.D0 - CALL self%partialDer(dPsi, dz, dr) + CALL self%partialDer(nNodes, dPsi, dz, dr) invJ(1,1:2) = (/ dr(2), -dz(2) /) invJ(2,1:2) = (/ -dr(1), dz(1) /) diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 3018989..705587a 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -40,10 +40,11 @@ MODULE moduleMesh3DCart END TYPE meshCell3DCart ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, dPsi, dx, dy, dz) + PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx, dy, dz) IMPORT meshCell3DCart CLASS(meshCell3DCart), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz END SUBROUTINE partialDer_interface @@ -59,7 +60,7 @@ MODULE moduleMesh3DCart !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() CONTAINS - PROCEDURE, PASS:: init => initCellTetra3DCart + PROCEDURE, PASS:: init => initCellTetra PROCEDURE, PASS:: randPos => randPosCellTetra PROCEDURE, PASS:: calcCell => volumeTetra PROCEDURE, PASS:: fPsi => fPsiTetra @@ -135,6 +136,7 @@ MODULE moduleMesh3DCart INTEGER:: s self%n = n + self%nNodes = SIZE(p) self%n1 => mesh%nodes(p(1))%obj self%n2 => mesh%nodes(p(2))%obj self%n3 => mesh%nodes(p(3))%obj @@ -170,13 +172,13 @@ MODULE moduleMesh3DCart END SUBROUTINE initEdge3DCartTria !Get nodes from surface - PURE FUNCTION getNodes3DCartTria(self) RESULT(n) + PURE FUNCTION getNodes3DCartTria(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshEdge3DCartTria), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - ALLOCATE(n(1:3)) n = (/self%n1%n, self%n2%n, self%n3%n/) END FUNCTION getNodes3DCartTria @@ -238,7 +240,7 @@ MODULE moduleMesh3DCart !VOLUME FUNCTIONS !TETRA FUNCTIONS !Inits tetrahedron element - SUBROUTINE initCellTetra3DCart(self, n, p, nodes) + SUBROUTINE initCellTetra(self, n, p, nodes) USE moduleRefParam IMPLICIT NONE @@ -282,7 +284,7 @@ MODULE moduleMesh3DCart ALLOCATE(self%listPart_in(1:nSpecies)) ALLOCATE(self%totalWeight(1:nSpecies)) - END SUBROUTINE initCellTetra3DCart + END SUBROUTINE initCellTetra !Random position in volume tetrahedron FUNCTION randPosCellTetra(self) RESULT(r) @@ -315,7 +317,7 @@ MODULE moduleMesh3DCart self%volume = 0.D0 Xi = (/0.25D0, 0.25D0, 0.25D0/) - self%volume = self%detJac(Xi) + self%volume = self%detJac(Xi, 4) END SUBROUTINE volumeTetra @@ -392,11 +394,12 @@ MODULE moduleMesh3DCart END FUNCTION dPsiTetraXi3 !Computes the derivatives in global coordinates - PURE SUBROUTINE partialDerTetra(self, dPsi, dx, dy, dz) + PURE SUBROUTINE partialDerTetra(self, nNodes, dPsi, dx, dy, dz) IMPLICIT NONE CLASS(meshCell3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: dPsi(1:3, 1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: dPsi(1:3, 1:nNodes) REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz dx(1) = DOT_PRODUCT(dPsi(1,:), self%x) @@ -413,11 +416,12 @@ MODULE moduleMesh3DCart END SUBROUTINE partialDerTetra - PURE FUNCTION elemKTetra(self) RESULT(localK) + PURE FUNCTION elemKTetra(self, nNodes) RESULT(localK) IMPLICIT NONE CLASS(meshCell3DCartTetra), INTENT(in):: self - REAL(8):: localK(1:self%nNodes,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) REAL(8):: invJ(1:3,1:3), detJ @@ -427,19 +431,20 @@ MODULE moduleMesh3DCart !TODO: One point Gauss integral. Upgrade when possible Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) dPsi = self%dPsi(Xi, 4) - detJ = self%detJac(Xi, dPsi) - invJ = self%invJac(Xi, dPsi) + detJ = self%detJac(Xi, 4, dPsi) + invJ = self%invJac(Xi, 4, dPsi) fPsi = self%fPsi(Xi, 4) localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ END FUNCTION elemKTetra - PURE FUNCTION elemFTetra(self, source) RESULT(localF) + PURE FUNCTION elemFTetra(self, nNodes, source) RESULT(localF) IMPLICIT NONE CLASS(meshCell3DCartTetra), INTENT(in):: self - REAL(8), INTENT(in):: source(1:self%nNodes) - REAL(8):: localF(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) REAL(8):: Xi(1:3) REAL(8):: detJ, f @@ -448,7 +453,7 @@ MODULE moduleMesh3DCart Xi = 0.D0 Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) dPsi = self%dPsi(Xi, 4) - detJ = self%detJac(Xi, dPsi) + detJ = self%detJac(Xi, 4, dPsi) fPsi = self%fPsi(Xi, 4) f = DOT_PRODUCT(fPsi, source) localF = f*fPsi*1.D0*detJ @@ -467,7 +472,7 @@ MODULE moduleMesh3DCart self%n3%emData%phi, & self%n4%emData%phi /) - array = -self%gatherDF(Xi, phi) + array = -self%gatherDF(Xi, 4, phi) END FUNCTION gatherEFTetra @@ -493,7 +498,7 @@ MODULE moduleMesh3DCart self%n3%emData%B(3), & self%n4%emData%B(3) /) - array = self%gatherF(Xi, 3, B) + array = self%gatherF(Xi, 4, B) END FUNCTION gatherMFTetra @@ -510,11 +515,12 @@ MODULE moduleMesh3DCart END FUNCTION insideTetra - PURE FUNCTION getNodesTetra(self) RESULT(n) + PURE FUNCTION getNodesTetra(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell3DCartTetra), INTENT(in):: self - INTEGER:: n(1:self%nnodes) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) @@ -533,8 +539,8 @@ MODULE moduleMesh3DCart Xi = 0.D0 deltaR = (/r(1) - self%x(1), r(2) - self%y(1), r(3) - self%z(1) /) dPsi = self%dPsi(Xi, 4) - invJ = self%invJac(Xi, dPsi) - detJ = self%detJac(Xi, dPsi) + invJ = self%invJac(Xi, 4, dPsi) + detJ = self%detJac(Xi, 4, dPsi) Xi = MATMUL(invJ, deltaR)/detJ END FUNCTION phy2logTetra @@ -567,14 +573,15 @@ 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, nNodes, dPsi_in) RESULT(dJ) IMPLICIT NONE CLASS(meshCell3DCart), INTENT(in)::self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3, 1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3, 1:nNodes) REAL(8):: dJ - REAL(8):: dPsi(1:3, 1:self%nNodes) + REAL(8):: dPsi(1:3, 1:nNodes) REAL(8):: dx(1:3), dy(1:3), dz(1:3) IF (PRESENT(dPsi_in)) THEN @@ -585,20 +592,21 @@ MODULE moduleMesh3DCart END IF - CALL self%partialDer(dPsi, dx, dy, dz) + CALL self%partialDer(nNodes, dPsi, dx, dy, dz) dJ = dx(1)*(dy(2)*dz(3) - dy(3)*dz(2)) & - dx(2)*(dy(1)*dz(3) - dy(3)*dz(1)) & + dx(3)*(dy(1)*dz(2) - dy(2)*dz(1)) END FUNCTION detJ3DCart - PURE FUNCTION invJ3DCart(self,Xi,dPsi_in) RESULT(invJ) + PURE FUNCTION invJ3DCart(self, Xi, nNodes, dPsi_in) RESULT(invJ) IMPLICIT NONE CLASS(meshCell3DCart), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3, 1:self%nNodes) - REAL(8):: dPsi(1:3, 1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3, 1:nNodes) + REAL(8):: dPsi(1:3, 1:nNodes) REAL(8), DIMENSION(1:3):: dx, dy, dz REAL(8):: invJ(1:3,1:3) @@ -610,7 +618,7 @@ MODULE moduleMesh3DCart END IF - CALL self%partialDer(dPsi, dx, dy, dz) + CALL self%partialDer(nNodes, dPsi, dx, dy, dz) invJ(1,1) = (dy(2)*dz(3) - dy(3)*dz(2)) invJ(1,2) = -(dy(1)*dz(3) - dy(3)*dz(1)) invJ(1,3) = (dy(1)*dz(2) - dy(2)*dz(1)) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index aeba4ba..6561850 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -66,6 +66,8 @@ MODULE moduleMesh !Parent of Edge element TYPE, PUBLIC, ABSTRACT, EXTENDS(meshElement):: meshEdge + !Nomber of nodes in the edge + INTEGER:: nNodes !Connectivity to cells CLASS(meshCell), POINTER:: e1 => NULL(), e2 => NULL() !Connectivity to cells in meshColl @@ -102,10 +104,11 @@ MODULE moduleMesh END SUBROUTINE initEdge_interface !Get nodes index from node - PURE FUNCTION getNodesEdge_interface(self) RESULT(n) + PURE FUNCTION getNodesEdge_interface(self, nNodes) RESULT(n) IMPORT:: meshEdge CLASS(meshEdge), INTENT(in):: self - INTEGER, ALLOCATABLE:: n(:) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) END FUNCTION getNodesEdge_interface @@ -166,7 +169,7 @@ MODULE moduleMesh !Init the cell PROCEDURE(initCell_interface), DEFERRED, PASS:: init !Get the index of the nodes - PROCEDURE(getNodesVol_interface), DEFERRED, PASS:: getNodes + PROCEDURE(getNodesCell_interface), DEFERRED, PASS:: getNodes !Calculate random position on the cell PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos !Obtain functions and values of cell natural functions @@ -208,12 +211,13 @@ MODULE moduleMesh END SUBROUTINE initCell_interface - PURE FUNCTION getNodesVol_interface(self) RESULT(n) + PURE FUNCTION getNodesCell_interface(self, nNodes) RESULT(n) IMPORT:: meshCell CLASS(meshCell), INTENT(in):: self - INTEGER:: n(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) - END FUNCTION getNodesVol_interface + END FUNCTION getNodesCell_interface PURE FUNCTION fPsi_interface(self, Xi, nNodes) RESULT(fPsi) IMPORT:: meshCell @@ -233,20 +237,22 @@ MODULE moduleMesh END FUNCTION dPsi_interface - PURE FUNCTION detJac_interface(self, Xi, dPsi_in) RESULT(dJ) + PURE FUNCTION detJac_interface(self, Xi, nNodes, dPsi_in) RESULT(dJ) IMPORT:: meshCell CLASS(meshCell), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) REAL(8):: dJ END FUNCTION detJac_interface - PURE FUNCTION invJac_interface(self, Xi, dPsi_in) RESULT(invJ) + PURE FUNCTION invJac_interface(self, Xi, nNodes, dPsi_in) RESULT(invJ) IMPORT:: meshCell CLASS(meshCell), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) REAL(8):: invJ(1:3,1:3) END FUNCTION invJac_interface @@ -259,18 +265,20 @@ MODULE moduleMesh END FUNCTION gatherArray_interface - PURE FUNCTION elemK_interface(self) RESULT(localK) + PURE FUNCTION elemK_interface(self, nNodes) RESULT(localK) IMPORT:: meshCell CLASS(meshCell), INTENT(in):: self - REAL(8):: localK(1:self%nNodes,1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) END FUNCTION elemK_interface - PURE FUNCTION elemF_interface(self, source) RESULT(localF) + PURE FUNCTION elemF_interface(self, nNodes, source) RESULT(localF) IMPORT:: meshCell CLASS(meshCell), INTENT(in):: self - REAL(8), INTENT(in):: source(1:self%nNodes) - REAL(8):: localF(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) END FUNCTION elemF_interface @@ -478,19 +486,22 @@ MODULE moduleMesh CONTAINS !Constructs the global K matrix - SUBROUTINE constructGlobalK(self) + PURE SUBROUTINE constructGlobalK(self) IMPLICIT NONE CLASS(meshParticles), INTENT(inout):: self INTEGER:: e + INTEGER:: nNodes INTEGER, ALLOCATABLE:: n(:) REAL(8), ALLOCATABLE:: localK(:,:) - INTEGER:: nNodes, i, j + INTEGER:: i, j DO e = 1, self%numCells - n = self%cells(e)%obj%getNodes() - localK = self%cells(e)%obj%elemK() - nNodes = SIZE(n) + nNodes = self%cells(e)%obj%nNodes + ALLOCATE(n(1:nNodes)) + ALLOCATE(localK(1:nNodes, 1:nNodes)) + n = self%cells(e)%obj%getNodes(nNodes) + localK = self%cells(e)%obj%elemK(nNodes) DO i = 1, nNodes DO j = 1, nNodes @@ -499,6 +510,8 @@ MODULE moduleMesh END DO END DO + + DEALLOCATE(n, localK) END DO @@ -523,51 +536,53 @@ MODULE moduleMesh END SUBROUTINE resetOutput !Gather the value of valNodes (scalar) at position Xi - PURE FUNCTION gatherF_scalar(self, Xi, valNodes) RESULT(f) + PURE FUNCTION gatherF_scalar(self, Xi, nNodes, valNodes) RESULT(f) IMPLICIT NONE CLASS(meshCell), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in):: valNodes(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: valNodes(1:nNodes) REAL(8):: f - REAL(8):: fPsi(1:self%nNodes) + REAL(8):: fPsi(1:nNodes) - fPsi = self%fPsi(Xi, self%nNodes) + fPsi = self%fPsi(Xi, nNodes) f = DOT_PRODUCT(fPsi, valNodes) END FUNCTION gatherF_scalar !Gather the value of valNodes (array) at position Xi - PURE FUNCTION gatherF_array(self, Xi, n, valNodes) RESULT(f) + PURE FUNCTION gatherF_array(self, Xi, nNodes, valNodes) RESULT(f) IMPLICIT NONE CLASS(meshCell), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: n - REAL(8), INTENT(in):: valNodes(1:self%nNodes, 1:n) - REAL(8):: f(1:n) - REAL(8):: fPsi(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: valNodes(1:nNodes, 1:3) + REAL(8):: f(1:3) + REAL(8):: fPsi(1:nNodes) - fPsi = self%fPsi(Xi, self%nNodes) + fPsi = self%fPsi(Xi, nNodes) f = MATMUL(fPsi, valNodes) END FUNCTION gatherF_array !Gather the spatial derivative of valNodes (scalar) at position Xi - PURE FUNCTION gatherDF_scalar(self, Xi, valNodes) RESULT(df) + PURE FUNCTION gatherDF_scalar(self, Xi, nNodes, valNodes) RESULT(df) IMPLICIT NONE CLASS(meshCell), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - REAL(8), INTENT(in):: valNodes(1:self%nNodes) + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: valNodes(1:nNodes) REAL(8):: df(1:3) - REAL(8):: dPsi(1:3, 1:self%nNodes) - REAL(8):: dPsiR(1:3, 1:self%nNodes) + REAL(8):: dPsi(1:3, 1:nNodes) + REAL(8):: dPsiR(1:3, 1:nNodes) REAL(8):: invJ(1:3, 1:3), detJ - dPsi = self%dPsi(Xi, self%nNodes) - detJ = self%detJac(Xi, dPsi) - invJ = self%invJac(Xi, dPsi) + dPsi = self%dPsi(Xi, nNodes) + detJ = self%detJac(Xi, nNodes, dPsi) + invJ = self%invJac(Xi, nNodes, dPsi) dPsiR = MATMUL(invJ, dPsi)/detJ df = (/ DOT_PRODUCT(dPsiR(1,:), valNodes), & DOT_PRODUCT(dPsiR(2,:), valNodes), & @@ -576,29 +591,30 @@ MODULE moduleMesh END FUNCTION gatherDF_scalar !Scatters particle properties into cell nodes - SUBROUTINE scatter(self, part) + SUBROUTINE scatter(self, nNodes, part) USE moduleMath USE moduleSpecies USE OMP_LIB IMPLICIT NONE CLASS(meshCell), INTENT(inout):: self + INTEGER, INTENT(in):: nNodes CLASS(particle), INTENT(in):: part - REAL(8):: fPsi(1:self%nNodes) - INTEGER:: cellNodes(1:self%nNodes) + REAL(8):: fPsi(1:nNodes) + INTEGER:: cellNodes(1:nNodes) REAL(8):: tensorS(1:3, 1:3) INTEGER:: sp INTEGER:: i CLASS(meshNode), POINTER:: node - cellNodes = self%getNodes() - fPsi = self%fPsi(part%Xi, self%nNodes) + cellNodes = self%getNodes(nNodes) + fPsi = self%fPsi(part%Xi, nNodes) tensorS = outerProduct(part%v, part%v) sp = part%species%n - DO i = 1, self%nNodes + DO i = 1, nNodes node => mesh%nodes(cellNodes(i))%obj CALL OMP_SET_LOCK(node%lock) node%output(sp)%den = node%output(sp)%den + part%weight*fPsi(i) diff --git a/src/modules/mesh/moduleMeshBoundary.f90 b/src/modules/mesh/moduleMeshBoundary.f90 index cc5d78c..713c091 100644 --- a/src/modules/mesh/moduleMeshBoundary.f90 +++ b/src/modules/mesh/moduleMeshBoundary.f90 @@ -55,10 +55,10 @@ MODULE moduleMeshBoundary !Scatter particle in associated volume IF (ASSOCIATED(edge%e1)) THEN - CALL edge%e1%scatter(part) + CALL edge%e1%scatter(edge%e1%nNodes, part) ELSE - CALL edge%e2%scatter(part) + CALL edge%e2%scatter(edge%e2%nNodes, part) END IF diff --git a/src/modules/solver/electromagnetic/moduleEM.f90 b/src/modules/solver/electromagnetic/moduleEM.f90 index 55eb618..bdf6b03 100644 --- a/src/modules/solver/electromagnetic/moduleEM.f90 +++ b/src/modules/solver/electromagnetic/moduleEM.f90 @@ -26,12 +26,14 @@ MODULE moduleEM CLASS(boundaryEM), INTENT(in):: self CLASS(meshEdge):: edge + INTEGER:: nNodes INTEGER, ALLOCATABLE:: nodes(:) INTEGER:: n - nodes = edge%getNodes() + nNodes = edge%nNodes + nodes = edge%getNodes(nNodes) - DO n = 1, SIZE(nodes) + DO n = 1, nNodes SELECT CASE(self%typeEM) CASE ("dirichlet") mesh%K(nodes(n), :) = 0.D0 @@ -66,8 +68,8 @@ MODULE moduleEM !$OMP DO REDUCTION(+:vectorF) DO e = 1, mesh%numCells - nodes = mesh%cells(e)%obj%getNodes() - nNodes = SIZE(nodes) + nNodes = mesh%cells(e)%obj%nNodes + nodes = mesh%cells(e)%obj%getNodes(nNodes) !Calculates charge density (rho) in element nodes ALLOCATE(rho(1:nNodes)) rho = 0.D0 @@ -79,7 +81,7 @@ MODULE moduleEM END DO !Calculates local F vector - localF = mesh%cells(e)%obj%elemF(rho) + localF = mesh%cells(e)%obj%elemF(nNodes, rho) !Assign local F to global F DO i = 1, nNodes diff --git a/src/modules/solver/moduleSolver.f90 b/src/modules/solver/moduleSolver.f90 index 6962075..02932c9 100644 --- a/src/modules/solver/moduleSolver.f90 +++ b/src/modules/solver/moduleSolver.f90 @@ -354,11 +354,13 @@ MODULE moduleSolver IMPLICIT NONE INTEGER:: n + CLASS(meshCell), POINTER:: cell !Loops over the particles to scatter them !$OMP DO DO n = 1, nPartOld - CALL mesh%cells(partOld(n)%vol)%obj%scatter(partOld(n)) + cell => mesh%cells(partOld(n)%vol)%obj + CALL cell%scatter(cell%nNodes, partOld(n)) END DO !$OMP END DO