diff --git a/src/modules/mesh/0D/moduleMesh0D.f90 b/src/modules/mesh/0D/moduleMesh0D.f90 index 5f89f20..133ae8e 100644 --- a/src/modules/mesh/0D/moduleMesh0D.f90 +++ b/src/modules/mesh/0D/moduleMesh0D.f90 @@ -6,7 +6,8 @@ MODULE moduleMesh0D TYPE, PUBLIC, EXTENDS(meshNode):: meshNode0D INTEGER:: n1 CONTAINS - PROCEDURE, PASS:: init => initNode0D + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode0D PROCEDURE, PASS:: getCoordinates => getCoord0D END TYPE meshNode0D @@ -14,20 +15,21 @@ MODULE moduleMesh0D TYPE, PUBLIC, EXTENDS(meshCell):: meshCell0D CLASS(meshNode), POINTER:: n1 CONTAINS - PROCEDURE, PASS:: init => initCell0D - PROCEDURE, PASS:: getNodes => getNodes0D - PROCEDURE, PASS:: randPos => randPos0D - PROCEDURE, PASS:: fPsi => fPsi0D - PROCEDURE, PASS:: dPsi => dPsi0D - PROCEDURE, PASS:: detJac => detJ0D - PROCEDURE, PASS:: invJac => invJ0D - PROCEDURE, PASS:: elemK => elemK0D - PROCEDURE, PASS:: elemF => elemF0D - PROCEDURE, PASS:: gatherElectricField => gatherEF0D - PROCEDURE, PASS:: gatherMagneticField => gatherMF0D - PROCEDURE, PASS:: phy2log => phy2log0D - PROCEDURE, NOPASS:: inside => inside0D - PROCEDURE, PASS:: nextElement => nextElement0D + PROCEDURE, PASS:: init => initCell0D + PROCEDURE, PASS:: getNodes => getNodes0D + PROCEDURE, PASS:: randPos => randPos0D + PROCEDURE, NOPASS:: fPsi => fPsi0D + PROCEDURE, NOPASS:: dPsi => dPsi0D + PROCEDURE, PASS:: partialDer => partialDer0D + PROCEDURE, NOPASS:: detJac => detJ0D + PROCEDURE, NOPASS:: invJac => invJ0D + PROCEDURE, PASS:: gatherElectricField => gatherEF0D + PROCEDURE, PASS:: gatherMagneticField => gatherMF0D + PROCEDURE, PASS:: elemK => elemK0D + PROCEDURE, PASS:: elemF => elemF0D + PROCEDURE, PASS:: phy2log => phy2log0D + PROCEDURE, NOPASS:: inside => inside0D + PROCEDURE, PASS:: neighbourElement => neighbourElement0D END TYPE meshCell0D @@ -89,6 +91,7 @@ MODULE moduleMesh0D END SUBROUTINE initCell0D + !Get the nodes of the volume PURE FUNCTION getNodes0D(self, nNodes) RESULT(n) IMPLICIT NONE @@ -100,6 +103,7 @@ MODULE moduleMesh0D END FUNCTION getNodes0D + !Calculate random position inside the volume FUNCTION randPos0D(self) RESULT(r) IMPLICIT NONE @@ -110,10 +114,9 @@ MODULE moduleMesh0D END FUNCTION randPos0D - PURE FUNCTION fPsi0D(self, Xi, nNodes) RESULT(fPsi) + PURE FUNCTION fPsi0D(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshCell0D), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) @@ -122,10 +125,9 @@ MODULE moduleMesh0D END FUNCTION fPsi0D - PURE FUNCTION dPsi0D(self, Xi, nNodes) RESULT(dPsi) + PURE FUNCTION dPsi0D(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - CLASS(meshCell0D), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: dPsi(1:3,1:nNodes) @@ -134,31 +136,17 @@ MODULE moduleMesh0D END FUNCTION dPsi0D - PURE FUNCTION detJ0D(self, Xi, nNodes, dPsi_in) RESULT(dJ) + PURE FUNCTION partialDer0D(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE CLASS(meshCell0D), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) - REAL(8):: dJ + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) + REAL(8):: pDer(1:3, 1:3) - dJ = 0.D0 + pDer = 0.D0 - END FUNCTION detJ0D - - PURE FUNCTION invJ0D(self, Xi, nNodes, dPsi_in) RESULT(invJ) - IMPLICIT NONE - - CLASS(meshCell0D), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - 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 + END FUNCTION partialDer0D PURE FUNCTION elemK0D(self, nNodes) RESULT(localK) IMPLICIT NONE @@ -234,15 +222,36 @@ MODULE moduleMesh0D END FUNCTION inside0D - SUBROUTINE nextElement0D(self, Xi, nextElement) + SUBROUTINE neighbourElement0D(self, Xi, neighbourElement) IMPLICIT NONE CLASS(meshCell0D), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement - nextElement => NULL() + neighbourElement => NULL() - END SUBROUTINE nextElement0D + END SUBROUTINE neighbourElement0D + + !COMMON FUNCTIONS + PURE FUNCTION detJ0D(pDer) RESULT(dJ) + IMPLICIT NONE + + REAL(8), INTENT(in):: pDer(1:3, 1:3) + REAL(8):: dJ + + dJ = 0.D0 + + END FUNCTION detJ0D + + PURE FUNCTION invJ0D(pDer) RESULT(invJ) + IMPLICIT NONE + + REAL(8), INTENT(in):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3) + + invJ = 0.D0 + + END FUNCTION invJ0D END MODULE moduleMesh0D diff --git a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 index 56d1b65..6a43ca1 100644 --- a/src/modules/mesh/1DCart/moduleMesh1DCart.f90 +++ b/src/modules/mesh/1DCart/moduleMesh1DCart.f90 @@ -14,7 +14,8 @@ MODULE moduleMesh1DCart !Element coordinates REAL(8):: x = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initNode1DCart + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode1DCart PROCEDURE, PASS:: getCoordinates => getCoord1DCart END TYPE meshNode1DCart @@ -25,6 +26,7 @@ MODULE moduleMesh1DCart !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL() CONTAINS + !meshEdge DEFERRED PROCEDURES PROCEDURE, PASS:: init => initEdge1DCart PROCEDURE, PASS:: getNodes => getNodes1DCart PROCEDURE, PASS:: intersection => intersection1DCart @@ -32,27 +34,7 @@ MODULE moduleMesh1DCart END TYPE meshEdge1DCart - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell1DCart - CONTAINS - PROCEDURE, PASS:: detJac => detJ1DCart - PROCEDURE, PASS:: invJac => invJ1DCart - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshCell1DCart - - ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx) - IMPORT meshCell1DCart - CLASS(meshCell1DCart), INTENT(in):: self - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) - REAL(8), INTENT(out), DIMENSION(1):: dx - - END SUBROUTINE partialDer_interface - - END INTERFACE - - TYPE, PUBLIC, EXTENDS(meshCell1DCart):: meshCell1DCartSegm + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell1DCartSegm !Element coordinates REAL(8):: x(1:2) !Connectivity to nodes @@ -61,20 +43,24 @@ MODULE moduleMesh1DCart CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() REAL(8):: arNodes(1:2) CONTAINS - PROCEDURE, PASS:: init => initCell1DCartSegm - PROCEDURE, PASS:: randPos => randPos1DCartSegm - PROCEDURE, PASS:: area => areaSegm - PROCEDURE, PASS:: fPsi => fPsiSegm - PROCEDURE, PASS:: dPsi => dPsiSegm - PROCEDURE, PASS:: partialDer => partialDerSegm - PROCEDURE, PASS:: elemK => elemKSegm - PROCEDURE, PASS:: elemF => elemFSegm - PROCEDURE, PASS:: gatherElectricField => gatherEFSegm - PROCEDURE, PASS:: gatherMagneticField => gatherMFSegm - PROCEDURE, NOPASS:: inside => insideSegm - PROCEDURE, PASS:: getNodes => getNodesSegm - PROCEDURE, PASS:: phy2log => phy2logSegm - PROCEDURE, PASS:: nextElement => nextElementSegm + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCell1DCartSegm + PROCEDURE, PASS:: getNodes => getNodesSegm + PROCEDURE, PASS:: randPos => randPos1DCartSegm + PROCEDURE, NOPASS:: fPsi => fPsiSegm + PROCEDURE, NOPASS:: dPsi => dPsiSegm + PROCEDURE, PASS:: partialDer => partialDerSegm + PROCEDURE, NOPASS:: detJac => detJ1DCart + PROCEDURE, NOPASS:: invJac => invJ1DCart + PROCEDURE, PASS:: gatherElectricField => gatherEFSegm + PROCEDURE, PASS:: gatherMagneticField => gatherMFSegm + PROCEDURE, PASS:: elemK => elemKSegm + PROCEDURE, PASS:: elemF => elemFSegm + PROCEDURE, NOPASS:: inside => insideSegm + PROCEDURE, PASS:: phy2log => phy2logSegm + PROCEDURE, PASS:: neighbourElement => neighbourElementSegm + !PARTICLUAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: area => areaSegm END TYPE meshCell1DCartSegm @@ -219,7 +205,19 @@ MODULE moduleMesh1DCart END SUBROUTINE initCell1DCartSegm - !Calculates a random position in 1D volume + !Get nodes from 1D volume + PURE FUNCTION getNodesSegm(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell1DCartSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) + + n = (/ self%n1%n, self%n2%n /) + + END FUNCTION getNodesSegm + + !Random position in 1D volume FUNCTION randPos1DCartSegm(self) RESULT(r) USE moduleRandom IMPLICIT NONE @@ -227,135 +225,63 @@ MODULE moduleMesh1DCart CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8):: r(1:3) REAL(8):: Xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) + REAL(8):: fPsi(1:2) - Xi(1) = random(-1.D0, 1.D0) - Xi(2:3) = 0.D0 + Xi = 0.D0 + Xi(1) = random(-1.D0, 1.D0) fPsi = self%fPsi(Xi, 2) + + r = 0.D0 r(1) = DOT_PRODUCT(fPsi, self%x) END FUNCTION randPos1DCartSegm - !Computes element area - PURE SUBROUTINE areaSegm(self) - IMPLICIT NONE - - CLASS(meshCell1DCartSegm), INTENT(inout):: self - REAL(8):: l !element length - REAL(8):: fPsi(1:2) - REAL(8):: detJ - REAL(8):: Xi(1:3) - - self%volume = 0.D0 - self%arNodes = 0.D0 - !1 point Gauss integral - Xi = 0.D0 - fPsi = self%fPsi(Xi, 2) - detJ = self%detJac(Xi, 2) - l = 2.D0*detJ - self%volume = l - self%arNodes = fPsi*l - - END SUBROUTINE areaSegm - !Computes element functions at point Xi - PURE FUNCTION fPsiSegm(self, Xi, nNodes) RESULT(fPsi) + PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) - fPsi(1) = 1.D0 - Xi(1) - fPsi(2) = 1.D0 + Xi(1) + fPsi = (/ 1.D0 - Xi(1), & + 1.D0 + Xi(1) /) - fPsi = fPsi * 5.D-1 + fPsi = fPsi * 0.50D0 END FUNCTION fPsiSegm - !Computes element derivative shape function at Xi - PURE FUNCTION dPsiSegm(self, Xi, nNodes) RESULT(dPsi) + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: dPsi(1:3,1:nNodes) dPsi = 0.D0 - dPsi(1, 1) = -5.D-1 - dPsi(1, 2) = 5.D-1 + dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /) END FUNCTION dPsiSegm - !Computes partial derivatives of coordinates - PURE SUBROUTINE partialDerSegm(self, nNodes, dPsi, dx) + !Partial derivative in global coordinates + PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self INTEGER, INTENT(in):: nNodes REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) - REAL(8), INTENT(out), DIMENSION(1):: dx + REAL(8):: pDer(1:3, 1:3) - dx(1) = DOT_PRODUCT(dPsi(1,:), self%x) + pDer = 0.D0 - END SUBROUTINE partialDerSegm + pDer(1,1) = DOT_PRODUCT(dPsi(1,1:2), self%x(1:2)) + pDer(2,2) = 1.D0 + pDer(3,3) = 1.D0 - !Computes local stiffness matrix - PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK) - IMPLICIT NONE - - CLASS(meshCell1DCartSegm), INTENT(in):: self - 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 - INTEGER:: l - - localK = 0.D0 - Xi = 0.D0 - DO l = 1, 3 - Xi(1) = corSeg(l) - dPsi = self%dPsi(Xi, 2) - 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 - - END DO - - END FUNCTION elemKSegm - - PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF) - IMPLICIT NONE - - CLASS(meshCell1DCartSegm), INTENT(in):: self - 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) - INTEGER:: l - - localF = 0.D0 - Xi = 0.D0 - - DO l = 1, 3 - Xi(1) = corSeg(l) - detJ = self%detJac(Xi, 2) - fPsi = self%fPsi(Xi, 2) - f = DOT_PRODUCT(fPsi, source) - localF = localF + f*fPsi*wSeg(l)*detJ - - END DO - - END FUNCTION elemFSegm + END FUNCTION partialDerSegm PURE FUNCTION gatherEFSegm(self, Xi) RESULT(array) IMPLICIT NONE @@ -391,6 +317,68 @@ MODULE moduleMesh1DCart END FUNCTION gatherMFSegm + !Computes element local stiffness matrix + PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK) + IMPLICIT NONE + + CLASS(meshCell1DCartSegm), INTENT(in):: self + 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):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3), detJ + INTEGER:: l + + localK = 0.D0 + + Xi = 0.D0 + !Start 1D Gauss Quad Integral + DO l = 1, 3 + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & + RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + wSeg(l)/detJ + + END DO + + END FUNCTION elemKSegm + + !Computes the local source vector for a force f + PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF) + IMPLICIT NONE + + CLASS(meshCell1DCartSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) + REAL(8):: fPsi(1:2) + REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3) + REAL(8):: Xi(1:3) + REAL(8):: detJ, f + INTEGER:: l + + localF = 0.D0 + + Xi = 0.D0 + !Start 1D Gauss Quad Integral + DO l = 1, 3 + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 2) + f = DOT_PRODUCT(fPsi, source) + localF = localF + f*fPsi*wSeg(l)*detJ + + END DO + + END FUNCTION elemFSegm + PURE FUNCTION insideSegm(Xi) RESULT(ins) IMPLICIT NONE @@ -402,101 +390,87 @@ MODULE moduleMesh1DCart END FUNCTION insideSegm - !Get nodes from 1D volume - PURE FUNCTION getNodesSegm(self, nNodes) RESULT(n) - IMPLICIT NONE - - CLASS(meshCell1DCartSegm), INTENT(in):: self - INTEGER, INTENT(in):: nNodes - INTEGER:: n(1:nNodes) - - n = (/ self%n1%n, self%n2%n /) - - END FUNCTION getNodesSegm - - PURE FUNCTION phy2logSegm(self, r) RESULT(xN) + PURE FUNCTION phy2logSegm(self, r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) + REAL(8):: Xi(1:3) - xN = 0.D0 - xN(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0 + Xi = 0.D0 + + Xi(1) = 2.D0*(r(1) - self%x(1))/(self%x(2) - self%x(1)) - 1.D0 END FUNCTION phy2logSegm - !Get next element for a logical position Xi - SUBROUTINE nextElementSegm(self, Xi, nextElement) + !Get the next element for a logical position Xi + SUBROUTINE neighbourElementSegm(self, Xi, neighbourElement) IMPLICIT NONE CLASS(meshCell1DCartSegm), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement - NULLIFY(nextElement) + NULLIFY(neighbourElement) IF (Xi(1) < -1.D0) THEN - nextElement => self%e2 + neighbourElement => self%e2 ELSEIF (Xi(1) > 1.D0) THEN - nextElement => self%e1 + neighbourElement => self%e1 END IF - END SUBROUTINE nextElementSegm + END SUBROUTINE neighbourElementSegm - !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS - !Calculates a random position in 1D volume - !Computes the element Jacobian determinant - PURE FUNCTION detJ1DCart(self, Xi, nNodes, dPsi_in) RESULT(dJ) + !Computes element area + PURE SUBROUTINE areaSegm(self) IMPLICIT NONE - CLASS(meshCell1DCart), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) - REAL(8):: dPsi(1:3,1:nNodes) + CLASS(meshCell1DCartSegm), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:2) + + self%volume = 0.D0 + self%arNodes = 0.D0 + !1D 1 point Gauss Quad Integral + Xi = 0.D0 + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 2) + !Computes total volume of the cell + self%volume = detJ*2.D0 + !Computes volume per node + self%arNodes = fPsi*self%volume + + END SUBROUTINE areaSegm + + !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS + !Computes element Jacobian determinant + PURE FUNCTION detJ1DCart(pDer) RESULT(dJ) + IMPLICIT NONE + + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: dJ - REAL(8):: dx(1) - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(Xi, 2) - - END IF - - CALL self%partialDer(2, dPsi, dx) - dJ = dx(1) + dJ = pDer(1, 1) END FUNCTION detJ1DCart - !Computes the invers Jacobian - PURE FUNCTION invJ1DCart(self, Xi, nNodes, dPsi_in) RESULT(invJ) + !Computes element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJ1DCart(pDer) RESULT(invJ) IMPLICIT NONE - CLASS(meshCell1DCart), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: invJ(1:3,1:3) - REAL(8):: dPsi(1:3,1:nNodes) - REAL(8):: dx(1) - - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(Xi, 2) - - END IF invJ = 0.D0 - CALL self%partialDer(2, dPsi, dx) - - invJ(1,1) = 1.D0/dx(1) + invJ(1, 1) = 1.D0/pDer(1, 1) + invJ(2, 2) = 1.D0 + invJ(3, 3) = 1.D0 END FUNCTION invJ1DCart diff --git a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 index 8b441be..8230901 100644 --- a/src/modules/mesh/1DRad/moduleMesh1DRad.f90 +++ b/src/modules/mesh/1DRad/moduleMesh1DRad.f90 @@ -8,13 +8,14 @@ MODULE moduleMesh1DRad IMPLICIT NONE REAL(8), PARAMETER:: corSeg(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /) - REAL(8), PARAMETER:: wSeg(1:3) = (/ 5.D0/9.D0 , 8.D0/9.D0, 5.D0/9.D0 /) + REAL(8), PARAMETER:: wSeg(1:3) = (/ 5.D0/9.D0 , 8.D0/9.D0, 5.D0/9.D0 /) TYPE, PUBLIC, EXTENDS(meshNode):: meshNode1DRad !Element coordinates REAL(8):: r = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initNode1DRad + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode1DRad PROCEDURE, PASS:: getCoordinates => getCoord1DRad END TYPE meshNode1DRad @@ -25,6 +26,7 @@ MODULE moduleMesh1DRad !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL() CONTAINS + !meshEdge DEFERRED PROCEDURES PROCEDURE, PASS:: init => initEdge1DRad PROCEDURE, PASS:: getNodes => getNodes1DRad PROCEDURE, PASS:: intersection => intersection1DRad @@ -32,28 +34,7 @@ MODULE moduleMesh1DRad END TYPE meshEdge1DRad - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell1DRad - CONTAINS - PROCEDURE, PASS:: detJac => detJ1DRad - PROCEDURE, PASS:: invJac => invJ1DRad - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshCell1DRad - - ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx) - IMPORT meshCell1DRad - - CLASS(meshCell1DRad), INTENT(in):: self - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) - REAL(8), INTENT(out), DIMENSION(1):: dx - - END SUBROUTINE partialDer_interface - - END INTERFACE - - TYPE, PUBLIC, EXTENDS(meshCell1DRad):: meshCell1DRadSegm + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell1DRadSegm !Element coordinates REAL(8):: r(1:2) !Connectivity to nodes @@ -62,20 +43,24 @@ MODULE moduleMesh1DRad CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL() REAL(8):: arNodes(1:2) CONTAINS - PROCEDURE, PASS:: init => initCell1DRadSegm - PROCEDURE, PASS:: randPos => randPos1DRadSeg - PROCEDURE, PASS:: area => areaRad - PROCEDURE, PASS:: fPsi => fPsiRad - PROCEDURE, PASS:: dPsi => dPsiRad - PROCEDURE, PASS:: partialDer => partialDerRad - PROCEDURE, PASS:: elemK => elemKRad - PROCEDURE, PASS:: elemF => elemFRad - PROCEDURE, PASS:: gatherElectricField => gatherEFRad - PROCEDURE, PASS:: gatherMagneticField => gatherMFRad - PROCEDURE, NOPASS:: inside => insideRad - PROCEDURE, PASS:: getNodes => getNodesRad - PROCEDURE, PASS:: phy2log => phy2logRad - PROCEDURE, PASS:: nextElement => nextElementRad + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCell1DRadSegm + PROCEDURE, PASS:: getNodes => getNodesSegm + PROCEDURE, PASS:: randPos => randPos1DRadSegm + PROCEDURE, NOPASS:: fPsi => fPsiSegm + PROCEDURE, NOPASS:: dPsi => dPsiSegm + PROCEDURE, PASS:: partialDer => partialDerSegm + PROCEDURE, NOPASS:: detJac => detJ1DRad + PROCEDURE, NOPASS:: invJac => invJ1DRad + PROCEDURE, PASS:: gatherElectricField => gatherEFSegm + PROCEDURE, PASS:: gatherMagneticField => gatherMFSegm + PROCEDURE, PASS:: elemK => elemKSegm + PROCEDURE, PASS:: elemF => elemFSegm + PROCEDURE, NOPASS:: inside => insideSegm + PROCEDURE, PASS:: phy2log => phy2logSegm + PROCEDURE, PASS:: neighbourElement => neighbourElementSegm + !PARTICLUAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: area => areaSegm END TYPE meshCell1DRadSegm @@ -139,7 +124,6 @@ MODULE moduleMesh1DRad self%r = r1(1) self%normal = (/ 1.D0, 0.D0, 0.D0 /) - self%normal = self%normal/NORM2(self%normal) !Boundary index self%boundary => boundary(bt) @@ -221,8 +205,20 @@ MODULE moduleMesh1DRad END SUBROUTINE initCell1DRadSegm - !Calculates a random position in 1D volume - FUNCTION randPos1DRadSeg(self) RESULT(r) + !Get nodes from 1D volume + PURE FUNCTION getNodesSegm(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell1DRadSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) + + n = (/ self%n1%n, self%n2%n /) + + END FUNCTION getNodesSegm + + !Random position in 1D volume + FUNCTION randPos1DRadSegm(self) RESULT(r) USE moduleRandom IMPLICIT NONE @@ -231,152 +227,63 @@ MODULE moduleMesh1DRad REAL(8):: Xi(1:3) REAL(8):: fPsi(1:2) - Xi(1) = random(-1.D0, 1.D0) - Xi(2:3) = 0.D0 + Xi = 0.D0 + Xi(1) = random(-1.D0, 1.D0) fPsi = self%fPsi(Xi, 2) + + r = 0.D0 r(1) = DOT_PRODUCT(fPsi, self%r) - END FUNCTION randPos1DRadSeg - - !Computes element area - PURE SUBROUTINE areaRad(self) - IMPLICIT NONE - - CLASS(meshCell1DRadSegm), INTENT(inout):: self - REAL(8):: l !element length - REAL(8):: fPsi(1:2), fPsi_node(1:2) - REAL(8):: r - REAL(8):: detJ - REAL(8):: Xi(1:3) - - self%volume = 0.D0 - self%arNodes = 0.D0 - !1 point Gauss integral - Xi = 0.D0 - fPsi = self%fPsi(Xi, 2) - 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, 2, self%r) - self%arNodes(1) = fPsi(1)*r*l - Xi = (/ 5.D-1, 0.D0, 0.D0/) - r = self%gatherF(Xi, 2, self%r) - self%arNodes(2) = fPsi(2)*r*l - - END SUBROUTINE areaRad + END FUNCTION randPos1DRadSegm !Computes element functions at point Xi - PURE FUNCTION fPsiRad(self, Xi, nNodes) RESULT(fPsi) + PURE FUNCTION fPsiSegm(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) - fPsi(1) = 1.D0 - Xi(1) - fPsi(2) = 1.D0 + Xi(1) + fPsi = (/ 1.D0 - Xi(1), & + 1.D0 + Xi(1) /) - fPsi = fPsi * 5.D-1 + fPsi = fPsi * 0.50D0 - END FUNCTION fPsiRad + END FUNCTION fPsiSegm - !Computes element derivative shape function at Xi - PURE FUNCTION dPsiRad(self, Xi, nNodes) RESULT(dPsi) + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiSegm(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: dPsi(1:3,1:nNodes) dPsi = 0.D0 - dPsi(1, 1) = -5.D-1 - dPsi(1, 2) = 5.D-1 + dPsi(1, 1:2) = (/ -5.D-1, 5.D-1 /) - END FUNCTION dPsiRad + END FUNCTION dPsiSegm - !Computes partial derivatives of coordinates - PURE SUBROUTINE partialDerRad(self, nNodes, dPsi, dx) + !Partial derivative in global coordinates + PURE FUNCTION partialDerSegm(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self INTEGER, INTENT(in):: nNodes REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) - REAL(8), INTENT(out), DIMENSION(1):: dx + REAL(8):: pDer(1:3, 1:3) - dx(1) = DOT_PRODUCT(dPsi(1,:), self%r) + pDer = 0.D0 - END SUBROUTINE partialDerRad + pDer(1,1) = DOT_PRODUCT(dPsi(1,1:2), self%r(1:2)) + pDer(2,2) = 1.D0 + pDer(3,3) = 1.D0 - !Computes local stiffness matrix - PURE FUNCTION elemKRad(self, nNodes) RESULT(localK) - USE moduleConstParam, ONLY: PI2 - IMPLICIT NONE + END FUNCTION partialDerSegm - CLASS(meshCell1DRadSegm), INTENT(in):: self - 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 - REAL(8):: r, fPsi(1:2) - INTEGER:: l - - localK = 0.D0 - Xi = 0.D0 - DO l = 1, 3 - Xi(1) = corSeg(l) - dPsi = self%dPsi(Xi, 2) - 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/)), & - RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & - r*wSeg(l)/detJ - - END DO - - localK = localK*PI2 - - END FUNCTION elemKRad - - PURE FUNCTION elemFRad(self, nNodes, source) RESULT(localF) - USE moduleConstParam, ONLY: PI2 - IMPLICIT NONE - - CLASS(meshCell1DRadSegm), INTENT(in):: self - 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) - INTEGER:: l - - localF = 0.D0 - Xi = 0.D0 - - DO l = 1, 3 - Xi(1) = corSeg(l) - detJ = self%detJac(Xi, 2) - fPsi = self%fPsi(Xi, 2) - r = DOT_PRODUCT(fPsi, self%r) - f = DOT_PRODUCT(fPsi, source) - localF = localF + f*fPsi*r*wSeg(l)*detJ - - END DO - - END FUNCTION elemFRad - - PURE FUNCTION gatherEFRad(self, Xi) RESULT(array) + PURE FUNCTION gatherEFSegm(self, Xi) RESULT(array) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) @@ -388,9 +295,9 @@ MODULE moduleMesh1DRad array = -self%gatherDF(Xi, 2, phi) - END FUNCTION gatherEFRad + END FUNCTION gatherEFSegm - PURE FUNCTION gatherMFRad(self, Xi) RESULT(array) + PURE FUNCTION gatherMFSegm(self, Xi) RESULT(array) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) @@ -408,9 +315,79 @@ MODULE moduleMesh1DRad array = self%gatherF(Xi, 2, B) - END FUNCTION gatherMFRad + END FUNCTION gatherMFSegm - PURE FUNCTION insideRad(Xi) RESULT(ins) + !Computes element local stiffness matrix + PURE FUNCTION elemKSegm(self, nNodes) RESULT(localK) + USE moduleConstParam, ONLY: PI2 + IMPLICIT NONE + + CLASS(meshCell1DRadSegm), INTENT(in):: self + 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):: pDer(1:3, 1:3) + REAL(8):: r + REAL(8):: invJ(1:3,1:3), detJ + INTEGER:: l + + localK = 0.D0 + + Xi = 0.D0 + !Start 1D Gauss Quad Integral + DO l = 1, 3 + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + r = self%gatherF(Xi, 4, self%r) + localK = localK + MATMUL(RESHAPE(MATMUL(invJ,dPsi), (/ 2, 1/)), & + RESHAPE(MATMUL(invJ,dPsi), (/ 1, 2/)))* & + r*wSeg(l)/detJ + + END DO + localK = localK*PI2 + + END FUNCTION elemKSegm + + !Computes the local source vector for a force f + PURE FUNCTION elemFSegm(self, nNodes, source) RESULT(localF) + USE moduleConstParam, ONLY: PI2 + IMPLICIT NONE + + CLASS(meshCell1DRadSegm), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) + REAL(8):: fPsi(1:2) + REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3) + REAL(8):: Xi(1:3) + REAL(8):: detJ, f + REAL(8):: r + INTEGER:: l + + localF = 0.D0 + + Xi = 0.D0 + !Start 1D Gauss Quad Integral + DO l = 1, 3 + Xi(1) = corSeg(l) + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 2) + r = DOT_PRODUCT(fPsi, self%r) + f = DOT_PRODUCT(fPsi, source) + localF = localF + r*f*fPsi*wSeg(l)*detJ + + END DO + localF = localF*PI2 + + END FUNCTION elemFSegm + + PURE FUNCTION insideSegm(Xi) RESULT(ins) IMPLICIT NONE REAL(8), INTENT(in):: Xi(1:3) @@ -419,102 +396,97 @@ MODULE moduleMesh1DRad ins = Xi(1) >=-1.D0 .AND. & Xi(1) <= 1.D0 - END FUNCTION insideRad + END FUNCTION insideSegm - !Get nodes from 1D volume - PURE FUNCTION getNodesRad(self, nNodes) RESULT(n) - IMPLICIT NONE - - CLASS(meshCell1DRadSegm), INTENT(in):: self - INTEGER, INTENT(in):: nNodes - INTEGER:: n(1:nNodes) - - n = (/ self%n1%n, self%n2%n /) - - END FUNCTION getNodesRad - - PURE FUNCTION phy2logRad(self, r) RESULT(xN) + PURE FUNCTION phy2logSegm(self, r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) - REAL(8):: xN(1:3) + REAL(8):: Xi(1:3) - xN = 0.D0 - xN(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.D0 + Xi = 0.D0 - END FUNCTION phy2logRad + Xi(1) = 2.D0*(r(1) - self%r(1))/(self%r(2) - self%r(1)) - 1.D0 - !Get next element for a logical position Xi - SUBROUTINE nextElementRad(self, Xi, nextElement) + END FUNCTION phy2logSegm + + !Get the next element for a logical position Xi + SUBROUTINE neighbourElementSegm(self, Xi, neighbourElement) IMPLICIT NONE CLASS(meshCell1DRadSegm), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement - NULLIFY(nextElement) + NULLIFY(neighbourElement) IF (Xi(1) < -1.D0) THEN - nextElement => self%e2 + neighbourElement => self%e2 ELSEIF (Xi(1) > 1.D0) THEN - nextElement => self%e1 + neighbourElement => self%e1 END IF - END SUBROUTINE nextElementRad + END SUBROUTINE neighbourElementSegm - !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS - !Computes the element Jacobian determinant - PURE FUNCTION detJ1DRad(self, Xi, nNodes, dPsi_in) RESULT(dJ) + !Computes element area + PURE SUBROUTINE areaSegm(self) + USE moduleConstParam, ONLY: PI IMPLICIT NONE - CLASS(meshCell1DRad), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) - REAL(8):: dPsi(1:3,1:nNodes) + CLASS(meshCell1DRadSegm), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:2), pDer(1:3, 1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:2) + REAL(8):: r + + self%volume = 0.D0 + self%arNodes = 0.D0 + !1D 1 point Gauss Quad Integral + Xi = 0.D0 + dPsi = self%dPsi(Xi, 2) + pDer = self%partialDer(2, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 2) + !Computes total volume of the cell + r = DOT_PRODUCT(fPsi, self%r) + self%volume = r*detJ*2.D0*PI !2PI + !Computes volume per node + Xi = (/-5.D-1, 0.D0, 0.D0/) + r = self%gatherF(Xi, 2, self%r) + self%arNodes(1) = fPsi(1)*self%volume + Xi = (/ 5.D-1, 0.D0, 0.D0/) + r = self%gatherF(Xi, 2, self%r) + self%arNodes(2) = fPsi(2)*self%volume + + END SUBROUTINE areaSegm + + !COMMON FUNCTIONS FOR 1D VOLUME ELEMENTS + !Computes element Jacobian determinant + PURE FUNCTION detJ1DRad(pDer) RESULT(dJ) + IMPLICIT NONE + + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: dJ - REAL(8):: dx(1) - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(Xi, 2) - - END IF - - CALL self%partialDer(nNodes, dPsi, dx) - dJ = dx(1) + dJ = pDer(1, 1) END FUNCTION detJ1DRad - !Computes the invers Jacobian - PURE FUNCTION invJ1DRad(self, Xi, nNodes, dPsi_in) RESULT(invJ) + !Computes element Jacobian inverse matrix (without determinant) + PURE FUNCTION invJ1DRad(pDer) RESULT(invJ) IMPLICIT NONE - CLASS(meshCell1DRad), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - 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), INTENT(in):: pDer(1:3, 1:3) REAL(8):: invJ(1:3,1:3) - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(Xi, 2) - - END IF - invJ = 0.D0 - CALL self%partialDer(nNodes, dPsi, dx) - - invJ(1,1) = 1.D0/dx(1) + invJ(1, 1) = 1.D0/pDer(1, 1) + invJ(2, 2) = 1.D0 + invJ(3, 3) = 1.D0 END FUNCTION invJ1DRad diff --git a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 index 2ccde7d..eab1266 100644 --- a/src/modules/mesh/2DCart/moduleMesh2DCart.f90 +++ b/src/modules/mesh/2DCart/moduleMesh2DCart.f90 @@ -19,7 +19,8 @@ MODULE moduleMesh2DCart !Element coordinates REAL(8):: x = 0.D0, y = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initNode2DCart + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode2DCart PROCEDURE, PASS:: getCoordinates => getCoord2DCart END TYPE meshNode2DCart @@ -30,35 +31,16 @@ MODULE moduleMesh2DCart !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CONTAINS - PROCEDURE, PASS:: init => initEdge2DCart - PROCEDURE, PASS:: getNodes => getNodes2DCart + !meshEdge DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initEdge2DCart + PROCEDURE, PASS:: getNodes => getNodes2DCart PROCEDURE, PASS:: intersection => intersection2DCartEdge - PROCEDURE, PASS:: randPos => randPosEdge + PROCEDURE, PASS:: randPos => randPosEdge END TYPE meshEdge2DCart - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell2DCart - CONTAINS - PROCEDURE, PASS:: detJac => detJ2DCart - PROCEDURE, PASS:: invJac => invJ2DCart - PROCEDURE(partialDer_interface), DEFERRED, PASS, PRIVATE:: partialDer - - END TYPE meshCell2DCart - - ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx, dy) - IMPORT meshCell2DCart - CLASS(meshCell2DCart), INTENT(in):: self - 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 - - END INTERFACE - !Quadrilateral volume element - TYPE, PUBLIC, EXTENDS(meshCell2DCart):: meshCell2DCartQuad + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCartQuad !Element coordinates REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0 !Connectivity to nodes @@ -68,25 +50,29 @@ MODULE moduleMesh2DCart REAL(8):: arNodes(1:4) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initCellQuad2DCart - PROCEDURE, PASS:: randPos => randPosCellQuad - PROCEDURE, PASS:: area => areaQuad - PROCEDURE, PASS:: fPsi => fPsiQuad - PROCEDURE, PASS:: dPsi => dPsiQuad - PROCEDURE, PASS, PRIVATE:: partialDer => partialDerQuad - PROCEDURE, PASS:: elemK => elemKQuad - PROCEDURE, PASS:: elemF => elemFQuad - PROCEDURE, PASS:: gatherElectricField => gatherEFQuad - PROCEDURE, PASS:: gatherMagneticField => gatherMFQuad - PROCEDURE, NOPASS:: inside => insideQuad - PROCEDURE, PASS:: getNodes => getNodesQuad - PROCEDURE, PASS:: phy2log => phy2logQuad - PROCEDURE, PASS:: nextElement => nextElementQuad + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCellQuad2DCart + PROCEDURE, PASS:: getNodes => getNodesQuad + PROCEDURE, PASS:: randPos => randPosCellQuad + PROCEDURE, NOPASS:: fPsi => fPsiQuad + PROCEDURE, NOPASS:: dPsi => dPsiQuad + PROCEDURE, PASS:: partialDer => partialDerQuad + PROCEDURE, NOPASS:: detJac => detJ2DCart + PROCEDURE, NOPASS:: invJac => invJ2DCart + PROCEDURE, PASS:: gatherElectricField => gatherEFQuad + PROCEDURE, PASS:: gatherMagneticField => gatherMFQuad + PROCEDURE, PASS:: elemK => elemKQuad + PROCEDURE, PASS:: elemF => elemFQuad + PROCEDURE, NOPASS:: inside => insideQuad + PROCEDURE, PASS:: phy2log => phy2logQuad + PROCEDURE, PASS:: neighbourElement => neighbourElementQuad + !PARTICLUAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: area => areaQuad END TYPE meshCell2DCartQuad !Triangular volume element - TYPE, PUBLIC, EXTENDS(meshCell2DCart):: meshCell2DCartTria + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCartTria !Element coordinates REAL(8):: x(1:3) = 0.D0, y(1:3) = 0.D0 !Connectivity to nodes @@ -96,20 +82,24 @@ MODULE moduleMesh2DCart REAL(8):: arNodes(1:3) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initCellTria2DCart - PROCEDURE, PASS:: randPos => randPosCellTria - PROCEDURE, PASS:: area => areaTria - PROCEDURE, PASS:: fPsi => fPsiTria - PROCEDURE, PASS:: dPsi => dPsiTria - PROCEDURE, PASS, PRIVATE:: partialDer => partialDerTria - PROCEDURE, PASS:: elemK => elemKTria - PROCEDURE, PASS:: elemF => elemFTria - PROCEDURE, PASS:: gatherElectricField => gatherEFTria - PROCEDURE, PASS:: gatherMagneticField => gatherMFTria - PROCEDURE, NOPASS:: inside => insideTria - PROCEDURE, PASS:: getNodes => getNodesTria - PROCEDURE, PASS:: phy2log => phy2logTria - PROCEDURE, PASS:: nextElement => nextElementTria + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCellTria2DCart + PROCEDURE, PASS:: getNodes => getNodesTria + PROCEDURE, PASS:: randPos => randPosCellTria + PROCEDURE, NOPASS:: fPsi => fPsiTria + PROCEDURE, NOPASS:: dPsi => dPsiTria + PROCEDURE, PASS:: partialDer => partialDerTria + PROCEDURE, NOPASS:: detJac => detJ2DCart + PROCEDURE, NOPASS:: invJac => invJ2DCart + PROCEDURE, PASS:: gatherElectricField => gatherEFTria + PROCEDURE, PASS:: gatherMagneticField => gatherMFTria + PROCEDURE, PASS:: elemK => elemKTria + PROCEDURE, PASS:: elemF => elemFTria + PROCEDURE, NOPASS:: inside => insideTria + PROCEDURE, PASS:: phy2log => phy2logTria + PROCEDURE, PASS:: neighbourElement => neighbourElementTria + !PARTICULAR PROCEDURES + PROCEDURE, PASS, PRIVATE:: area => areaTria END TYPE meshCell2DCartTria @@ -175,6 +165,7 @@ MODULE moduleMesh2DCart r2 = self%n2%getCoordinates() self%x = (/r1(1), r2(1)/) self%y = (/r1(2), r2(2)/) + self%weight = 1.D0 !Normal vector self%normal = (/ -(self%y(2)-self%y(1)), & self%x(2)-self%x(1) , & @@ -290,84 +281,17 @@ MODULE moduleMesh2DCart END SUBROUTINE initCellQuad2DCart - !Computes element area - PURE SUBROUTINE areaQuad(self) - IMPLICIT NONE - - CLASS(meshCell2DCartQuad), INTENT(inout):: self - 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)*4.D0 !4 - fPsi = self%fPsi(Xi, 4) - self%volume = detJ - self%arNodes = fPsi*detJ - - END SUBROUTINE areaQuad - - !Computes element functions in point Xi - PURE FUNCTION fPsiQuad(self, Xi, nNodes) RESULT(fPsi) - IMPLICIT NONE - - CLASS(meshCell2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8):: fPsi(1:nNodes) - - fPsi(1) = (1.D0-Xi(1)) * (1.D0-Xi(2)) - fPsi(2) = (1.D0+Xi(1)) * (1.D0-Xi(2)) - fPsi(3) = (1.D0+Xi(1)) * (1.D0+Xi(2)) - fPsi(4) = (1.D0-Xi(1)) * (1.D0+Xi(2)) - - fPsi = fPsi*0.25D0 - - END FUNCTION fPsiQuad - - !Derivative element function at coordinates Xi - PURE FUNCTION dPsiQuad(self, Xi, nNodes) RESULT(dPsi) - IMPLICIT NONE - - CLASS(meshCell2DCartQuad), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8):: dPsi(1:3,1:nNodes) - - dPsi = 0.D0 - - dPsi(1,:) = (/ -(1.D0 - Xi(2)), & - (1.D0 - Xi(2)), & - (1.D0 + Xi(2)), & - -(1.D0 + Xi(2)) /) - - dPsi(2,:) = (/ -(1.D0 - Xi(1)), & - -(1.D0 + Xi(1)), & - (1.D0 + Xi(1)), & - (1.D0 - Xi(1)) /) - - dPsi = dPsi * 0.25D0 - - END FUNCTION dPsiQuad - - !Partial derivative in global coordinates - PURE SUBROUTINE partialDerQuad(self, nNodes, dPsi, dx, dy) + !Get nodes from quadrilateral element + PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) - REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy + INTEGER:: n(1:nNodes) - dx = (/ DOT_PRODUCT(dPsi(1,1:4),self%x(1:4)), & - DOT_PRODUCT(dPsi(2,1:4),self%x(1:4)) /) - dy = (/ DOT_PRODUCT(dPsi(1,1:4),self%y(1:4)), & - DOT_PRODUCT(dPsi(2,1:4),self%y(1:4)) /) + n = (/ self%n1%n, self%n2%n, self%n3%n, self%n4%n /) - END SUBROUTINE partialDerQuad + END FUNCTION getNodesQuad !Random position in quadrilateral volume FUNCTION randPosCellQuad(self) RESULT(r) @@ -379,78 +303,77 @@ MODULE moduleMesh2DCart REAL(8):: Xi(1:3) REAL(8):: fPsi(1:4) + Xi = 0.D0 Xi(1) = random(-1.D0, 1.D0) Xi(2) = random(-1.D0, 1.D0) - Xi(3) = 0.D0 fPsi = self%fPsi(Xi, 4) + r = 0.D0 r(1) = DOT_PRODUCT(fPsi, self%x) r(2) = DOT_PRODUCT(fPsi, self%y) - r(3) = 0.D0 END FUNCTION randPosCellQuad - !Computes element local stiffness matrix - PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK) + !Computes element functions in point Xi + PURE FUNCTION fPsiQuad(Xi, nNodes) RESULT(fPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: fPsi(1:nNodes) + + fPsi = (/ (1.D0-Xi(1)) * (1.D0-Xi(2)), & + (1.D0+Xi(1)) * (1.D0-Xi(2)), & + (1.D0+Xi(1)) * (1.D0+Xi(2)), & + (1.D0-Xi(1)) * (1.D0+Xi(2)) /) + + fPsi = fPsi * 0.25D0 + + END FUNCTION fPsiQuad + + !Derivative element function at coordinates Xi + PURE FUNCTION dPsiQuad(Xi, nNodes) RESULT(dPsi) + IMPLICIT NONE + + REAL(8), INTENT(in):: Xi(1:3) + INTEGER, INTENT(in):: nNodes + REAL(8):: dPsi(1:3,1:nNodes) + + dPsi = 0.D0 + + dPsi(1, 1:4) = (/ -(1.D0 - Xi(2)), & + (1.D0 - Xi(2)), & + (1.D0 + Xi(2)), & + -(1.D0 + Xi(2)) /) + + dPsi(2, 1:4) = (/ -(1.D0 - Xi(1)), & + -(1.D0 + Xi(1)), & + (1.D0 + Xi(1)), & + (1.D0 - Xi(1)) /) + + dPsi = dPsi * 0.25D0 + + END FUNCTION dPsiQuad + + !Partial derivative in global coordinates + PURE FUNCTION partialDerQuad(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self 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 - INTEGER:: l, m + REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) + REAL(8):: pDer(1:3, 1:3) - localK=0.D0 - Xi=0.D0 - !Start 2D Gauss Quad Integral - DO l=1, 3 - Xi(2) = corQuad(l) - DO m = 1, 3 - Xi(1) = corQuad(m) - fPsi = self%fPsi(Xi, 4) - dPsi = self%dPsi(Xi, 4) - 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 + pDer = 0.D0 - END DO - END DO + pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%x(1:4)), & + DOT_PRODUCT(dPsi(2,1:4),self%x(1:4)) /) + pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%y(1:4)), & + DOT_PRODUCT(dPsi(2,1:4),self%y(1:4)) /) + pDer(3,3) = 1.D0 - END FUNCTION elemKQuad - - !Computes the local source vector for a force f - PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF) - IMPLICIT NONE - - CLASS(meshCell2DCartQuad), INTENT(in):: self - 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 - INTEGER:: l, m - - localF = 0.D0 - Xi = 0.D0 - DO l=1, 3 - Xi(1) = corQuad(l) - DO m = 1, 3 - Xi(2) = corQuad(m) - detJ = self%detJac(Xi, 4) - fPsi = self%fPsi(Xi, 4) - f = DOT_PRODUCT(fPsi,source) - localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ - - END DO - END DO - - END FUNCTION elemFQuad + END FUNCTION partialDerQuad PURE FUNCTION gatherEFQuad(self, Xi) RESULT(array) IMPLICIT NONE @@ -494,6 +417,75 @@ MODULE moduleMesh2DCart END FUNCTION gatherMFQuad + !Computes element local stiffness matrix + PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK) + IMPLICIT NONE + + CLASS(meshCell2DCartQuad), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: r + REAL(8):: invJ(1:3,1:3), detJ + INTEGER:: l, m + + localK = 0.D0 + + Xi = 0.D0 + !Start 2D Gauss Quad Integral + DO l = 1, 3 + Xi(2) = corQuad(l) + DO m = 1, 3 + Xi(1) = corQuad(m) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + MATMUL(invJ,dPsi))* & + wQuad(l)*wQuad(m)/detJ + + END DO + END DO + + END FUNCTION elemKQuad + + !Computes the local source vector for a force f + PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF) + IMPLICIT NONE + + CLASS(meshCell2DCartQuad), INTENT(in):: self + 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):: dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: detJ, f + INTEGER:: l, m + + localF = 0.D0 + + Xi = 0.D0 + DO l = 1, 3 + Xi(1) = corQuad(l) + DO m = 1, 3 + Xi(2) = corQuad(m) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 4) + f = DOT_PRODUCT(fPsi,source) + localF = localF + f*fPsi*wQuad(l)*wQuad(m)*detJ + + END DO + END DO + + END FUNCTION elemFQuad + !Checks if a particle is inside a quad element PURE FUNCTION insideQuad(Xi) RESULT(ins) IMPLICIT NONE @@ -506,18 +498,6 @@ MODULE moduleMesh2DCart END FUNCTION insideQuad - !Gets nodes from quadrilateral element - PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n) - IMPLICIT NONE - - CLASS(meshCell2DCartQuad), INTENT(in):: self - INTEGER, INTENT(in):: nNodes - INTEGER:: n(1:nNodes) - - n = (/self%n1%n, self%n2%n, self%n3%n, self%n4%n /) - - END FUNCTION getNodesQuad - !Transforms physical coordinates to element coordinates PURE FUNCTION phy2logQuad(self,r) RESULT(Xi) IMPLICIT NONE @@ -527,6 +507,7 @@ MODULE moduleMesh2DCart REAL(8):: Xi(1:3) REAL(8):: XiO(1:3), detJ, invJ(1:3,1:3), f(1:3) REAL(8):: dPsi(1:3,1:4), fPsi(1:4) + REAL(8):: pDer(1:3, 1:3) REAL(8):: conv !Iterative newton method to transform coordinates @@ -535,7 +516,9 @@ MODULE moduleMesh2DCart DO WHILE(conv > 1.D-2) dPsi = self%dPsi(XiO, 4) - invJ = self%invJac(XiO, 4, dPsi) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) fPsi = self%fPsi(XiO, 4) f = (/ DOT_PRODUCT(fPsi,self%x), & DOT_PRODUCT(fPsi,self%y), & @@ -550,31 +533,56 @@ MODULE moduleMesh2DCart END FUNCTION phy2logQuad !Gets the next element for a logical position Xi - SUBROUTINE nextElementQuad(self, Xi, nextElement) + SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement) IMPLICIT NONE CLASS(meshCell2DCartQuad), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement REAL(8):: XiArray(1:4) INTEGER:: nextInt 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) + NULLIFY(neighbourElement) SELECT CASE (nextInt) CASE (1) - nextElement => self%e1 + neighbourElement => self%e1 CASE (2) - nextElement => self%e2 + neighbourElement => self%e2 CASE (3) - nextElement => self%e3 + neighbourElement => self%e3 CASE (4) - nextElement => self%e4 + neighbourElement => self%e4 END SELECT - END SUBROUTINE nextElementQuad + END SUBROUTINE neighbourElementQuad + + !Computes element area + PURE SUBROUTINE areaQuad(self) + IMPLICIT NONE + + CLASS(meshCell2DCartQuad), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: detJ + REAL(8):: fPsi(1:4) + REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) + + self%volume = 0.D0 + self%arNodes = 0.D0 + !2D 1 point Gauss Quad Integral + Xi = 0.D0 + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 4) + !Computes total volume of the cell + self%volume = detJ + !Computes volume per node + self%arNodes = fPsi*detJ + + END SUBROUTINE areaQuad !TRIA ELEMENT !Init tria element @@ -617,6 +625,18 @@ MODULE moduleMesh2DCart END SUBROUTINE initCellTria2DCart + !Gets node indexes from triangular element + PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell2DCartTria), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + INTEGER:: n(1:nNodes) + + n = (/self%n1%n, self%n2%n, self%n3%n /) + + END FUNCTION getNodesTria + !Random position in quadrilateral volume FUNCTION randPosCellTria(self) RESULT(r) USE moduleRandom @@ -639,31 +659,10 @@ MODULE moduleMesh2DCart END FUNCTION randPosCellTria - !Calculates area for triangular element - PURE SUBROUTINE areaTria(self) - IMPLICIT NONE - - CLASS(meshCell2DCartTria), INTENT(inout):: self - 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, 4)/2.D0 - fPsi = self%fPsi(Xi, 4) - self%volume = detJ - self%arNodes = fPsi*detJ - - END SUBROUTINE areaTria - !Shape functions for triangular element - PURE FUNCTION fPsiTria(self, Xi, nNodes) RESULT(fPsi) + PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) @@ -675,10 +674,9 @@ MODULE moduleMesh2DCart END FUNCTION fPsiTria !Derivative element function at coordinates Xi - PURE FUNCTION dPsiTria(self, Xi, nNodes) RESULT(dPsi) + PURE FUNCTION dPsiTria(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: dPsi(1:3,1:nNodes) @@ -690,76 +688,22 @@ MODULE moduleMesh2DCart END FUNCTION dPsiTria - PURE SUBROUTINE partialDerTria(self, nNodes, dPsi, dx, dy) + PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE CLASS(meshCell2DCartTria), INTENT(in):: self INTEGER, INTENT(in):: nNodes REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) - REAL(8), INTENT(out), DIMENSION(1:2):: dx, dy + REAL(8):: pDer(1:3, 1:3) - dx = (/ DOT_PRODUCT(dPsi(1,:),self%x), & - DOT_PRODUCT(dPsi(2,:),self%x) /) - dy = (/ DOT_PRODUCT(dPsi(1,:),self%y), & - DOT_PRODUCT(dPsi(2,:),self%y) /) + pDer = 0.D0 - END SUBROUTINE partialDerTria + pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%x(1:3)), & + DOT_PRODUCT(dPsi(2,1:3),self%x(1:3)) /) + pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%y(1:3)), & + DOT_PRODUCT(dPsi(2,1:3),self%y(1:3)) /) - !Computes element local stiffness matrix - PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) - IMPLICIT NONE - - CLASS(meshCell2DCartTria), INTENT(in):: self - 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 - INTEGER:: l - - localK=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, 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 - - END FUNCTION elemKTria - - !Computes element local source vector - PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF) - IMPLICIT NONE - - CLASS(meshCell2DCartTria), INTENT(in):: self - 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 - INTEGER:: l - - localF = 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, 3) - fPsi = self%fPsi(Xi, 3) - f = DOT_PRODUCT(fPsi,source) - localF = localF + f*fPsi*wTria(l)*detJ - - END DO - - END FUNCTION elemFTria + END FUNCTION partialDerTria PURE FUNCTION gatherEFTria(self, Xi) RESULT(array) IMPLICIT NONE @@ -799,6 +743,66 @@ MODULE moduleMesh2DCart END FUNCTION gatherMFTria + !Computes element local stiffness matrix + PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) + IMPLICIT NONE + + CLASS(meshCell2DCartTria), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8):: localK(1:nNodes,1:nNodes) + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3,1:3) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3), detJ + INTEGER:: l + + localK=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, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*wTria(l)/detJ + + END DO + + END FUNCTION elemKTria + + !Computes element local source vector + PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF) + IMPLICIT NONE + + CLASS(meshCell2DCartTria), INTENT(in):: self + INTEGER, INTENT(in):: nNodes + REAL(8), INTENT(in):: source(1:nNodes) + REAL(8):: localF(1:nNodes) + REAL(8):: fPsi(1:3) + REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3) + REAL(8):: Xi(1:3) + REAL(8):: detJ, f + INTEGER:: l + + localF = 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, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 3) + f = DOT_PRODUCT(fPsi,source) + localF = localF + f*fPsi*wTria(l)*detJ + + END DO + + END FUNCTION elemFTria + PURE FUNCTION insideTria(Xi) RESULT(ins) IMPLICIT NONE @@ -811,18 +815,6 @@ MODULE moduleMesh2DCart END FUNCTION insideTria - !Gets node indexes from triangular element - PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) - IMPLICIT NONE - - CLASS(meshCell2DCartTria), INTENT(in):: self - INTEGER, INTENT(in):: nNodes - INTEGER:: n(1:nNodes) - - n = (/self%n1%n, self%n2%n, self%n3%n /) - - END FUNCTION getNodesTria - !Transforms physical coordinates to element coordinates PURE FUNCTION phy2logTria(self,r) RESULT(Xi) IMPLICIT NONE @@ -830,96 +822,94 @@ MODULE moduleMesh2DCart CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: r(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:3) + REAL(8):: dPsi(1:3, 1:3) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3, 1:3), detJ !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, 3) - invJ = self%invJac(Xi, 3, dPsi) - detJ = self%detJac(Xi, 3, dPsi) + pDer = self%partialDer(3, dPsi) + invJ = self%invJac(pDer) + detJ = self%detJac(pDer) Xi = MATMUL(invJ,deltaR)/detJ END FUNCTION phy2logTria - SUBROUTINE nextElementTria(self, Xi, nextElement) + SUBROUTINE neighbourElementTria(self, Xi, neighbourElement) IMPLICIT NONE CLASS(meshCell2DCartTria), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement REAL(8):: XiArray(1:3) INTEGER:: nextInt XiArray = (/ Xi(2), 1.D0-Xi(1)-Xi(2), Xi(1) /) nextInt = MINLOC(XiArray,1) - NULLIFY(nextElement) + NULLIFY(neighbourElement) SELECT CASE (nextInt) CASE (1) - nextElement => self%e1 + neighbourElement => self%e1 CASE (2) - nextElement => self%e2 + neighbourElement => self%e2 CASE (3) - nextElement => self%e3 + neighbourElement => self%e3 END SELECT - END SUBROUTINE nextElementTria + END SUBROUTINE neighbourElementTria + + !Calculates area for triangular element + PURE SUBROUTINE areaTria(self) + IMPLICIT NONE + + CLASS(meshCell2DCartTria), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 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 /) + dPsi = self%dPsi(Xi, 3) + pDer = self%partialDer(3, dPsi) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 4) + !Computes total volume of the cell + self%volume = detJ + !Computes volume per node + self%arNodes = fPsi*detJ + + END SUBROUTINE areaTria !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 2D !Computes element Jacobian determinant - PURE FUNCTION detJ2DCart(self, Xi, nNodes, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ2DCart(pDer) RESULT(dJ) IMPLICIT NONE - CLASS(meshCell2DCart), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: dJ - REAL(8):: dPsi(1:3,1:nNodes) - REAL(8):: dx(1:2), dy(1:2) - IF(PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(Xi, 4) - - END IF - - CALL self%partialDer(nNodes, dPsi, dx, dy) - - dJ = dx(1)*dy(2)-dx(2)*dy(1) + dJ = pDer(1,1)*pDer(2,2)-pDer(1,2)*pDer(2,1) END FUNCTION detJ2DCart !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJ2DCart(self, Xi, nNodes, dPsi_in) RESULT(invJ) + PURE FUNCTION invJ2DCart(pDer) RESULT(invJ) IMPLICIT NONE - CLASS(meshCell2DCart), INTENT(in):: self - REAL(8), INTENT(in):: Xi(1:3) - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: invJ(1:3,1:3) - REAL(8):: dPsi(1:3,1:nNodes) - REAL(8):: dx(1:2), dy(1:2) - - IF(PRESENT(dPsi_in)) THEN - dPsi=dPsi_in - - ELSE - dPsi = self%dPsi(Xi, 4) - - END IF invJ = 0.D0 - CALL self%partialDer(nNodes, dPsi, dx, dy) - - invJ(1,1:2) = (/ dy(2), -dx(2) /) - invJ(2,1:2) = (/ -dy(1), dx(1) /) + invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) + invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(3, 3) = 1.D0 END FUNCTION invJ2DCart diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index 1e8c79d..d9a7f32 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -67,7 +67,7 @@ MODULE moduleMesh2DCyl PROCEDURE, PASS:: phy2log => phy2logQuad PROCEDURE, PASS:: neighbourElement => neighbourElementQuad !PARTICLUAR PROCEDURES - PROCEDURE, PASS:: area => areaQuad + PROCEDURE, PASS, PRIVATE:: area => areaQuad END TYPE meshCell2DCylQuad @@ -99,7 +99,7 @@ MODULE moduleMesh2DCyl PROCEDURE, PASS:: phy2log => phy2logTria PROCEDURE, PASS:: neighbourElement => neighbourElementTria !PARTICULAR PROCEDURES - PROCEDURE, PASS:: area => areaTria + PROCEDURE, PASS, PRIVATE:: area => areaTria END TYPE meshCell2DCylTria @@ -256,8 +256,13 @@ MODULE moduleMesh2DCyl 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 @@ -428,7 +433,7 @@ MODULE moduleMesh2DCyl 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):: dPsi(1:3, 1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: r REAL(8):: invJ(1:3,1:3), detJ @@ -445,8 +450,7 @@ MODULE moduleMesh2DCyl pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) invJ = self%invJac(pDer) - fPsi = self%fPsi(Xi, 4) - r = DOT_PRODUCT(fPsi,self%r) + r = self%gatherF(Xi, 4, self%r) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & MATMUL(invJ,dPsi))* & r*wQuad(l)*wQuad(m)/detJ @@ -467,7 +471,8 @@ MODULE moduleMesh2DCyl REAL(8), INTENT(in):: source(1:nNodes) REAL(8):: localF(1:nNodes) REAL(8):: Xi(1:3) - REAL(8):: fPsi(1:4), dPsi(1:3, 1:4) + REAL(8):: fPsi(1:4) + REAL(8):: dPsi(1:3, 1:4) REAL(8):: pDer(1:3, 1:3) REAL(8):: r REAL(8):: detJ, f @@ -483,7 +488,7 @@ MODULE moduleMesh2DCyl pDer = self%partialDer(4, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 4) - r = DOT_PRODUCT(fPsi,self%r) + r = DOT_PRODUCT(fPsi, self%r) f = DOT_PRODUCT(fPsi,source) localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ @@ -539,7 +544,7 @@ MODULE moduleMesh2DCyl END FUNCTION phy2logQuad - !Gets the next element for a logical position Xi + !Get the next element for a logical position Xi SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement) IMPLICIT NONE @@ -572,7 +577,8 @@ MODULE moduleMesh2DCyl IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(inout):: self - REAL(8):: r, Xi(1:3) + REAL(8):: Xi(1:3) + REAL(8):: r REAL(8):: detJ REAL(8):: fPsi(1:4) REAL(8):: dPsi(1:3, 1:4), pDer(1:3, 1:3) @@ -583,24 +589,24 @@ MODULE moduleMesh2DCyl Xi = 0.D0 dPsi = self%dPsi(Xi, 4) pDer = self%partialDer(4, dPsi) - detJ = self%detJac(pDer)*PI8 !4*2*pi + detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 4) !Computes total volume of the cell r = DOT_PRODUCT(fPsi,self%r) - self%volume = r*detJ + self%volume = r*detJ*PI8 !4*2*pi !Computes volume per node Xi = (/-5.D-1, -5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) - self%arNodes(1) = fPsi(1)*r*detJ + self%arNodes(1) = fPsi(1)*self%volume Xi = (/ 5.D-1, -5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) - self%arNodes(2) = fPsi(2)*r*detJ + self%arNodes(2) = fPsi(2)*self%volume Xi = (/ 5.D-1, 5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) - self%arNodes(3) = fPsi(3)*r*detJ + self%arNodes(3) = fPsi(3)*self%volume Xi = (/-5.D-1, 5.D-1, 0.D0/) r = self%gatherF(Xi, 4, self%r) - self%arNodes(4) = fPsi(4)*r*detJ + self%arNodes(4) = fPsi(4)*self%volume END SUBROUTINE areaQuad @@ -619,7 +625,7 @@ MODULE moduleMesh2DCyl !Assign node index self%n = n - !Assign nomber of nodes to cell + !Assign number of nodes of cell self%nNodes = SIZE(p) !Assign nodes to element @@ -736,7 +742,7 @@ MODULE moduleMesh2DCyl self%n2%emData%phi, & self%n3%emData%phi /) - array = -self%gatherDF(Xi, 4, phi) + array = -self%gatherDF(Xi, 3, phi) END FUNCTION gatherEFTria @@ -772,8 +778,8 @@ MODULE moduleMesh2DCyl INTEGER, INTENT(in):: nNodes REAL(8):: localK(1:nNodes,1:nNodes) REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3,1:3) REAL(8):: r - REAL(8):: fPsi(1:3), dPsi(1:3,1:3) REAL(8):: pDer(1:3, 1:3) REAL(8):: invJ(1:3,1:3), detJ INTEGER:: l @@ -788,8 +794,7 @@ MODULE moduleMesh2DCyl pDer = self%partialDer(3, dPsi) detJ = self%detJac(pDer) invJ = self%invJac(pDer) - fPsi = self%fPsi(Xi, 3) - r = DOT_PRODUCT(fPsi,self%r) + r = self%gatherF(Xi, 3, self%r) localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ END DO @@ -809,22 +814,23 @@ MODULE moduleMesh2DCyl REAL(8):: fPsi(1:3) REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3) REAL(8):: Xi(1:3) - REAL(8):: r REAL(8):: detJ, f + REAL(8):: r INTEGER:: l localF = 0.D0 - Xi = 0.D0 + + Xi = 0.D0 !Start 2D Gauss Quad Integral - DO l=1, 4 + DO l = 1, 4 Xi(1) = Xi1Tria(l) Xi(2) = Xi2Tria(l) dPsi = self%dPsi(Xi, 3) pDer = self%partialDer(3, dPsi) detJ = self%detJac(pDer) fPsi = self%fPsi(Xi, 3) - r = DOT_PRODUCT(fPsi,self%r) - f = DOT_PRODUCT(fPsi,source) + r = DOT_PRODUCT(fPsi, self%r) + f = DOT_PRODUCT(fPsi, source) localF = localF + r*f*fPsi*wTria(l)*detJ END DO @@ -897,10 +903,10 @@ MODULE moduleMesh2DCyl CLASS(meshCell2DCylTria), INTENT(inout):: self REAL(8):: Xi(1:3) - REAL(8):: r REAL(8):: dPsi(1:3, 1:3), pDer(1:3, 1:3) REAL(8):: detJ REAL(8):: fPsi(1:3) + REAL(8):: r self%volume = 0.D0 self%arNodes = 0.D0 @@ -908,13 +914,13 @@ MODULE moduleMesh2DCyl Xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /) dPsi = self%dPsi(Xi, 3) pDer = self%partialDer(3, dPsi) - detJ = self%detJac(pDer)*PI !2PI*1/2 - fPsi = self%fPsi(Xi, 4) + detJ = self%detJac(pDer) + fPsi = self%fPsi(Xi, 3) !Computes total volume of the cell - r = DOT_PRODUCT(fPsi,self%r) - self%volume = r*detJ + r = DOT_PRODUCT(fPsi, self%r) + self%volume = r*detJ*PI !2PI*1/2 !Computes volume per node - self%arNodes = fPsi*r*detJ + self%arNodes = fPsi*self%volume END SUBROUTINE areaTria @@ -939,9 +945,9 @@ MODULE moduleMesh2DCyl invJ = 0.D0 - invJ(1,1:2) = (/ pDer(2,2), -pDer(1,2) /) - invJ(2,1:2) = (/ -pDer(2,1), pDer(1,1) /) - invJ(3,3) = 1.D0 + invJ(1, 1:2) = (/ pDer(2,2), -pDer(1,2) /) + invJ(2, 1:2) = (/ -pDer(2,1), pDer(1,1) /) + invJ(3, 3) = 1.D0 END FUNCTION invJ2DCyl diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index a2c849d..34474cd 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -30,7 +30,7 @@ MODULE moduleMesh3DCart PROCEDURE, PASS:: intersection => intersection3DCartTria PROCEDURE, PASS:: randPos => randPosEdgeTria !PARTICULAR PROCEDURES - PROCEDURE, NOPASS:: fPsi => fPsiEdgeTria + PROCEDURE, NOPASS, PRIVATE:: fPsi => fPsiEdgeTria END TYPE meshEdge3DCartTria @@ -60,7 +60,7 @@ MODULE moduleMesh3DCart PROCEDURE, PASS:: phy2log => phy2logTetra PROCEDURE, PASS:: neighbourElement => neighbourElementTetra !PARTICULAR PROCEDURES - PROCEDURE, PASS:: calcVol => volumeTetra + PROCEDURE, PASS, PRIVATE:: calcVol => volumeTetra END TYPE meshCell3DCartTetra