From ba272de4e3c711e0a13706283ebcf1c6ae96a464 Mon Sep 17 00:00:00 2001 From: JGonzalez Date: Fri, 6 Jan 2023 12:16:54 +0100 Subject: [PATCH] DOES NOT COMPILE: Break Small break of changing functions. Still some geometries to change. --- src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 | 747 +++++++++---------- src/modules/mesh/3DCart/moduleMesh3DCart.f90 | 370 ++++----- src/modules/mesh/moduleMesh.f90 | 152 ++-- 3 files changed, 589 insertions(+), 680 deletions(-) diff --git a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 index d33fbbf..1e8c79d 100644 --- a/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 +++ b/src/modules/mesh/2DCyl/moduleMesh2DCyl.f90 @@ -19,7 +19,8 @@ MODULE moduleMesh2DCyl !Element coordinates REAL(8):: r = 0.D0, z = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initNode2DCyl + !meshNode DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initNode2DCyl PROCEDURE, PASS:: getCoordinates => getCoord2DCyl END TYPE meshNode2DCyl @@ -30,35 +31,16 @@ MODULE moduleMesh2DCyl !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL() CONTAINS - PROCEDURE, PASS:: init => initEdge2DCyl - PROCEDURE, PASS:: getNodes => getNodes2DCyl + !meshEdge DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initEdge2DCyl + PROCEDURE, PASS:: getNodes => getNodes2DCyl PROCEDURE, PASS:: intersection => intersection2DCylEdge - PROCEDURE, PASS:: randPos => randPosEdge + PROCEDURE, PASS:: randPos => randPosEdge END TYPE meshEdge2DCyl - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell2DCyl - CONTAINS - PROCEDURE, PASS:: detJac => detJ2DCyl - PROCEDURE, PASS:: invJac => invJ2DCyl - PROCEDURE(partialDer_interface), DEFERRED, PASS, PRIVATE:: partialDer - - END TYPE meshCell2DCyl - - ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dz, dr) - IMPORT meshCell2DCyl - CLASS(meshCell2DCyl), INTENT(in):: self - 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 - - END INTERFACE - !Quadrilateral volume element - TYPE, PUBLIC, EXTENDS(meshCell2DCyl):: meshCell2DCylQuad + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylQuad !Element coordinates REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0 !Connectivity to nodes @@ -68,25 +50,29 @@ MODULE moduleMesh2DCyl REAL(8):: arNodes(1:4) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initCellQuad2DCyl - 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 => initCellQuad2DCyl + PROCEDURE, PASS:: getNodes => getNodesQuad + PROCEDURE, PASS:: randPos => randPosCellQuad + PROCEDURE, NOPASS:: fPsi => fPsiQuad + PROCEDURE, NOPASS:: dPsi => dPsiQuad + PROCEDURE, PASS:: partialDer => partialDerQuad + PROCEDURE, NOPASS:: detJac => detJ2DCyl + PROCEDURE, NOPASS:: invJac => invJ2DCyl + 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:: area => areaQuad END TYPE meshCell2DCylQuad !Triangular volume element - TYPE, PUBLIC, EXTENDS(meshCell2DCyl):: meshCell2DCylTria + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylTria !Element coordinates REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0 !Connectivity to nodes @@ -96,20 +82,24 @@ MODULE moduleMesh2DCyl REAL(8):: arNodes(1:3) = 0.D0 CONTAINS - PROCEDURE, PASS:: init => initCellTria2DCyl - 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 => initCellTria2DCyl + PROCEDURE, PASS:: getNodes => getNodesTria + PROCEDURE, PASS:: randPos => randPosCellTria + PROCEDURE, NOPASS:: fPsi => fPsiTria + PROCEDURE, NOPASS:: dPsi => dPsiTria + PROCEDURE, PASS:: partialDer => partialDerTria + PROCEDURE, NOPASS:: detJac => detJ2DCyl + PROCEDURE, NOPASS:: invJac => invJ2DCyl + 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:: area => areaTria END TYPE meshCell2DCylTria @@ -294,99 +284,17 @@ MODULE moduleMesh2DCyl END SUBROUTINE initCellQuad2DCyl - !Computes element area - PURE SUBROUTINE areaQuad(self) - USE moduleConstParam, ONLY: PI8 - IMPLICIT NONE - - CLASS(meshCell2DCylQuad), INTENT(inout):: self - REAL(8):: r, Xi(1:3) - REAL(8):: detJ - REAL(8):: fPsi(1:4), fPsi_node(1:4) - - self%volume = 0.D0 - self%arNodes = 0.D0 - !2D 1 point Gauss Quad Integral - Xi = 0.D0 - detJ = self%detJac(Xi, 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, 4, self%r) - self%arNodes(1) = fPsi(1)*r*detJ - Xi = (/ 5.D-1, -5.D-1, 0.D0/) - 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, 4, self%r) - self%arNodes(3) = fPsi(3)*r*detJ - Xi = (/-5.D-1, 5.D-1, 0.D0/) - r = self%gatherF(Xi, 4, self%r) - self%arNodes(4) = fPsi(4)*r*detJ - - END SUBROUTINE areaQuad - - !Computes element functions in point Xi - PURE FUNCTION fPsiQuad(self, Xi, nNodes) RESULT(fPsi) - IMPLICIT NONE - - CLASS(meshCell2DCylQuad), INTENT(in):: self - 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(self, Xi, nNodes) RESULT(dPsi) - IMPLICIT NONE - - CLASS(meshCell2DCylQuad), 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, dz, dr) + !Gets nodes from quadrilateral element + PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n) IMPLICIT NONE CLASS(meshCell2DCylQuad), INTENT(in):: self INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) - REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr + INTEGER:: n(1:nNodes) - dz = (/ DOT_PRODUCT(dPsi(1,1:4),self%z(1:4)), & - DOT_PRODUCT(dPsi(2,1:4),self%z(1:4)) /) - dr = (/ DOT_PRODUCT(dPsi(1,1:4),self%r(1:4)), & - DOT_PRODUCT(dPsi(2,1:4),self%r(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) @@ -410,74 +318,64 @@ MODULE moduleMesh2DCyl END FUNCTION randPosCellQuad - !Computes element local stiffness matrix - PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK) - USE moduleConstParam, ONLY: PI2 + !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.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 FUNCTION partialDerQuad(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE CLASS(meshCell2DCylQuad), 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):: r - 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) - r = DOT_PRODUCT(fPsi,self%r) - localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & - MATMUL(invJ,dPsi))* & - r*wQuad(l)*wQuad(m)/detJ + pDer = 0.D0 - END DO - END DO - localK = localK*PI2 + pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%z(1:4)), & + DOT_PRODUCT(dPsi(2,1:4),self%z(1:4)) /) + pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:4),self%r(1:4)), & + DOT_PRODUCT(dPsi(2,1:4),self%r(1:4)) /) - END FUNCTION elemKQuad - - !Computes the local source vector for a force f - PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF) - USE moduleConstParam, ONLY: PI2 - IMPLICIT NONE - - CLASS(meshCell2DCylQuad), 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):: r - 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) - r = DOT_PRODUCT(fPsi,self%r) - f = DOT_PRODUCT(fPsi,source) - localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ - - END DO - END DO - localF = localF*PI2 - - END FUNCTION elemFQuad + END FUNCTION partialDerQuad PURE FUNCTION gatherEFQuad(self, Xi) RESULT(array) IMPLICIT NONE @@ -521,6 +419,80 @@ MODULE moduleMesh2DCyl END FUNCTION gatherMFQuad + !Computes element local stiffness matrix + PURE FUNCTION elemKQuad(self, nNodes) RESULT(localK) + USE moduleConstParam, ONLY: PI2 + IMPLICIT NONE + + CLASS(meshCell2DCylQuad), 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):: 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) + fPsi = self%fPsi(Xi, 4) + r = DOT_PRODUCT(fPsi,self%r) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)), & + MATMUL(invJ,dPsi))* & + r*wQuad(l)*wQuad(m)/detJ + + END DO + END DO + localK = localK*PI2 + + END FUNCTION elemKQuad + + !Computes the local source vector for a force f + PURE FUNCTION elemFQuad(self, nNodes, source) RESULT(localF) + USE moduleConstParam, ONLY: PI2 + IMPLICIT NONE + + CLASS(meshCell2DCylQuad), 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), dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: r + 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) + r = DOT_PRODUCT(fPsi,self%r) + f = DOT_PRODUCT(fPsi,source) + localF = localF + r*f*fPsi*wQuad(l)*wQuad(m)*detJ + + END DO + END DO + localF = localF*PI2 + + END FUNCTION elemFQuad + !Checks if a particle is inside a quad element PURE FUNCTION insideQuad(Xi) RESULT(ins) IMPLICIT NONE @@ -533,18 +505,6 @@ MODULE moduleMesh2DCyl END FUNCTION insideQuad - !Gets nodes from quadrilateral element - PURE FUNCTION getNodesQuad(self, nNodes) RESULT(n) - IMPLICIT NONE - - CLASS(meshCell2DCylQuad), 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 @@ -554,6 +514,7 @@ MODULE moduleMesh2DCyl 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 @@ -562,8 +523,9 @@ MODULE moduleMesh2DCyl DO WHILE(conv > 1.D-2) dPsi = self%dPsi(XiO, 4) - invJ = self%invJac(XiO, 4, dPsi) - detJ = self%detJac(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%z), & DOT_PRODUCT(fPsi,self%r), & @@ -578,31 +540,69 @@ MODULE moduleMesh2DCyl 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(meshCell2DCylQuad), 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) + USE moduleConstParam, ONLY: PI8 + IMPLICIT NONE + + CLASS(meshCell2DCylQuad), INTENT(inout):: self + REAL(8):: r, 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)*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, 4, self%r) + self%arNodes(1) = fPsi(1)*r*detJ + Xi = (/ 5.D-1, -5.D-1, 0.D0/) + 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, 4, self%r) + self%arNodes(3) = fPsi(3)*r*detJ + Xi = (/-5.D-1, 5.D-1, 0.D0/) + r = self%gatherF(Xi, 4, self%r) + self%arNodes(4) = fPsi(4)*r*detJ + + END SUBROUTINE areaQuad !TRIA ELEMENT !Init tria element @@ -645,6 +645,18 @@ MODULE moduleMesh2DCyl END SUBROUTINE initCellTria2DCyl + !Gets node indexes from triangular element + PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell2DCylTria), 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 @@ -667,36 +679,10 @@ MODULE moduleMesh2DCyl END FUNCTION randPosCellTria - !Calculates area for triangular element - PURE SUBROUTINE areaTria(self) - USE moduleConstParam, ONLY: PI - IMPLICIT NONE - - CLASS(meshCell2DCylTria), INTENT(inout):: self - REAL(8):: Xi(1:3) - REAL(8):: r - 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, 3)*PI !2PI*1/2 - 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 - self%arNodes = fPsi*r*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(meshCell2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) @@ -708,10 +694,9 @@ MODULE moduleMesh2DCyl 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(meshCell2DCylTria), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: dPsi(1:3,1:nNodes) @@ -723,84 +708,22 @@ MODULE moduleMesh2DCyl END FUNCTION dPsiTria - PURE SUBROUTINE partialDerTria(self, nNodes, dPsi, dz, dr) + PURE FUNCTION partialDerTria(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE CLASS(meshCell2DCylTria), INTENT(in):: self INTEGER, INTENT(in):: nNodes REAL(8), INTENT(in):: dPsi(1:3,1:nNodes) - REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr + REAL(8):: pDer(1:3, 1:3) - dz = (/ DOT_PRODUCT(dPsi(1,:),self%z), & - DOT_PRODUCT(dPsi(2,:),self%z) /) - dr = (/ DOT_PRODUCT(dPsi(1,:),self%r), & - DOT_PRODUCT(dPsi(2,:),self%r) /) + pDer = 0.D0 - END SUBROUTINE partialDerTria + pDer(1, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%z(1:3)), & + DOT_PRODUCT(dPsi(2,1:3),self%z(1:3)) /) + pDer(2, 1:2) = (/ DOT_PRODUCT(dPsi(1,1:3),self%r(1:3)), & + DOT_PRODUCT(dPsi(2,1:3),self%r(1:3)) /) - !Computes element local stiffness matrix - PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) - USE moduleConstParam, ONLY: PI2 - IMPLICIT NONE - - CLASS(meshCell2DCylTria), INTENT(in):: self - 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) - 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, 4) - 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 - - END DO - localK = localK*PI2 - - END FUNCTION elemKTria - - !Computes element local source vector - PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF) - USE moduleConstParam, ONLY: PI2 - IMPLICIT NONE - - CLASS(meshCell2DCylTria), 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):: r - 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) - r = DOT_PRODUCT(fPsi,self%r) - f = DOT_PRODUCT(fPsi,source) - localF = localF + r*f*fPsi*wTria(l)*detJ - - END DO - localF = localF*PI2 - - END FUNCTION elemFTria + END FUNCTION partialDerTria PURE FUNCTION gatherEFTria(self, Xi) RESULT(array) IMPLICIT NONE @@ -840,6 +763,75 @@ MODULE moduleMesh2DCyl END FUNCTION gatherMFTria + !Computes element local stiffness matrix + PURE FUNCTION elemKTria(self, nNodes) RESULT(localK) + USE moduleConstParam, ONLY: PI2 + IMPLICIT NONE + + CLASS(meshCell2DCylTria), INTENT(in):: self + 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) + 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) + fPsi = self%fPsi(Xi, 3) + r = DOT_PRODUCT(fPsi,self%r) + localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ + + END DO + localK = localK*PI2 + + END FUNCTION elemKTria + + !Computes element local source vector + PURE FUNCTION elemFTria(self, nNodes, source) RESULT(localF) + USE moduleConstParam, ONLY: PI2 + IMPLICIT NONE + + CLASS(meshCell2DCylTria), 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):: r + 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) + r = DOT_PRODUCT(fPsi,self%r) + f = DOT_PRODUCT(fPsi,source) + localF = localF + r*f*fPsi*wTria(l)*detJ + + END DO + localF = localF*PI2 + + END FUNCTION elemFTria + PURE FUNCTION insideTria(Xi) RESULT(ins) IMPLICIT NONE @@ -852,18 +844,6 @@ MODULE moduleMesh2DCyl END FUNCTION insideTria - !Gets node indexes from triangular element - PURE FUNCTION getNodesTria(self, nNodes) RESULT(n) - IMPLICIT NONE - - CLASS(meshCell2DCylTria), 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 @@ -871,96 +851,97 @@ MODULE moduleMesh2DCyl CLASS(meshCell2DCylTria), 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%z(1), r(2) - self%r(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(meshCell2DCylTria), 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) + USE moduleConstParam, ONLY: PI + IMPLICIT NONE + + 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) + + 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)*PI !2PI*1/2 + 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 + self%arNodes = fPsi*r*detJ + + END SUBROUTINE areaTria !COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS !Computes element Jacobian determinant - PURE FUNCTION detJ2DCyl(self, Xi, nNodes, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ2DCyl(pDer) RESULT(dJ) IMPLICIT NONE - CLASS(meshCell2DCyl), 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):: dz(1:2), dr(1:2) - IF(PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(Xi, nNodes) - - END IF - - CALL self%partialDer(nNodes, dPsi, dz, dr) - - dJ = dz(1)*dr(2)-dz(2)*dr(1) + dJ = pDer(1,1)*pDer(2,2)-pDer(1,2)*pDer(2,1) END FUNCTION detJ2DCyl !Computes element Jacobian inverse matrix (without determinant) - PURE FUNCTION invJ2DCyl(self, Xi, nNodes, dPsi_in) RESULT(invJ) + PURE FUNCTION invJ2DCyl(pDer) RESULT(invJ) IMPLICIT NONE - CLASS(meshCell2DCyl), 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):: dz(1:2), dr(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, dz, dr) - - invJ(1,1:2) = (/ dr(2), -dz(2) /) - invJ(2,1:2) = (/ -dr(1), dz(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 invJ2DCyl diff --git a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 index 705587a..a2c849d 100644 --- a/src/modules/mesh/3DCart/moduleMesh3DCart.f90 +++ b/src/modules/mesh/3DCart/moduleMesh3DCart.f90 @@ -11,6 +11,7 @@ MODULE moduleMesh3DCart !Element coordinates REAL(8):: x, y, z CONTAINS + !meshNode DEFERRED PROCEDURES PROCEDURE, PASS:: init => initNode3DCart PROCEDURE, PASS:: getCoordinates => getCoord3DCart @@ -23,36 +24,18 @@ MODULE moduleMesh3DCart !Connectivity to nodes CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL() CONTAINS + !meshEdge DEFERRED PROCEDURES PROCEDURE, PASS:: init => initEdge3DCartTria PROCEDURE, PASS:: getNodes => getNodes3DCartTria PROCEDURE, PASS:: intersection => intersection3DCartTria PROCEDURE, PASS:: randPos => randPosEdgeTria + !PARTICULAR PROCEDURES PROCEDURE, NOPASS:: fPsi => fPsiEdgeTria END TYPE meshEdge3DCartTria - TYPE, PUBLIC, ABSTRACT, EXTENDS(meshCell):: meshCell3DCart - CONTAINS - PROCEDURE, PASS:: detJac => detJ3DCart - PROCEDURE, PASS:: invJac => invJ3DCart - PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer - - END TYPE meshCell3DCart - - ABSTRACT INTERFACE - PURE SUBROUTINE partialDer_interface(self, nNodes, dPsi, dx, dy, dz) - IMPORT meshCell3DCart - CLASS(meshCell3DCart), INTENT(in):: self - 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 - - END INTERFACE - !Tetrahedron volume element - TYPE, PUBLIC, EXTENDS(meshCell3DCart):: meshCell3DCartTetra + TYPE, PUBLIC, EXTENDS(meshCell):: meshCell3DCartTetra !Element Coordinates REAL(8):: x(1:4) = 0.D0, y(1:4) = 0.D0, z(1:4) = 0.D0 !Connectivity to nodes @@ -60,22 +43,24 @@ MODULE moduleMesh3DCart !Connectivity to adjacent elements CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL() CONTAINS - PROCEDURE, PASS:: init => initCellTetra - PROCEDURE, PASS:: randPos => randPosCellTetra - PROCEDURE, PASS:: calcCell => volumeTetra - PROCEDURE, PASS:: fPsi => fPsiTetra - PROCEDURE, PASS:: dPsi => dPsiTetra - PROCEDURE, NOPASS, PRIVATE:: dPsiXi1 => dPsiTetraXi1 - PROCEDURE, NOPASS, PRIVATE:: dPsiXi2 => dPsiTetraXi2 - PROCEDURE, PASS:: partialDer => partialDerTetra - PROCEDURE, PASS:: elemK => elemKTetra - PROCEDURE, PASS:: elemF => elemFTetra - PROCEDURE, PASS:: gatherElectricField => gatherEFTetra - PROCEDURE, PASS:: gatherMagneticField => gatherMFTetra - PROCEDURE, NOPASS:: inside => insideTetra - PROCEDURE, PASS:: getNodes => getNodesTetra - PROCEDURE, PASS:: phy2log => phy2logTetra - PROCEDURE, PASS:: nextElement => nextElementTetra + !meshCell DEFERRED PROCEDURES + PROCEDURE, PASS:: init => initCellTetra + PROCEDURE, PASS:: getNodes => getNodesTetra + PROCEDURE, PASS:: randPos => randPosCellTetra + PROCEDURE, NOPASS:: fPsi => fPsiTetra + PROCEDURE, NOPASS:: dPsi => dPsiTetra + PROCEDURE, PASS:: partialDer => partialDerTetra + PROCEDURE, NOPASS:: detJac => detJ3DCart + PROCEDURE, NOPASS:: invJac => invJ3DCart + PROCEDURE, PASS:: gatherElectricField => gatherEFTetra + PROCEDURE, PASS:: gatherMagneticField => gatherMFTetra + PROCEDURE, PASS:: elemK => elemKTetra + PROCEDURE, PASS:: elemF => elemFTetra + PROCEDURE, NOPASS:: inside => insideTetra + PROCEDURE, PASS:: phy2log => phy2logTetra + PROCEDURE, PASS:: neighbourElement => neighbourElementTetra + !PARTICULAR PROCEDURES + PROCEDURE, PASS:: calcVol => volumeTetra END TYPE meshCell3DCartTetra @@ -227,13 +212,11 @@ MODULE moduleMesh3DCart IMPLICIT NONE REAL(8), INTENT(in):: Xi(1:3) - REAL(8), ALLOCATABLE:: fPsi(:) - - ALLOCATE(fPsi(1:3)) + REAL(8):: fPsi(1:3) fPsi(1) = 1.D0 - Xi(1) - Xi(2) fPsi(2) = Xi(1) - fPsi(3) = Xi(2) + fPsi(3) = Xi(2) END FUNCTION fPsiEdgeTria @@ -268,7 +251,7 @@ MODULE moduleMesh3DCart self%z = (/r1(3), r2(3), r3(3), r4(3)/) !Computes the element volume - CALL self%calcCell() + CALL self%calcVol() !Assign proportional volume to each node Xi = (/0.25D0, 0.25D0, 0.25D0/) @@ -286,6 +269,17 @@ MODULE moduleMesh3DCart END SUBROUTINE initCellTetra + PURE FUNCTION getNodesTetra(self, nNodes) RESULT(n) + IMPLICIT NONE + + CLASS(meshCell3DCartTetra), 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 getNodesTetra + !Random position in volume tetrahedron FUNCTION randPosCellTetra(self) RESULT(r) USE moduleRandom @@ -308,24 +302,10 @@ MODULE moduleMesh3DCart END FUNCTION randPosCellTetra - !Computes the element volume - PURE SUBROUTINE volumeTetra(self) - IMPLICIT NONE - - CLASS(meshCell3DCartTetra), INTENT(inout):: self - REAL(8):: Xi(1:3) - - self%volume = 0.D0 - Xi = (/0.25D0, 0.25D0, 0.25D0/) - self%volume = self%detJac(Xi, 4) - - END SUBROUTINE volumeTetra - !Computes element functions in point Xi - PURE FUNCTION fPsiTetra(self, Xi, nNodes) RESULT(fPsi) + PURE FUNCTION fPsiTetra(Xi, nNodes) RESULT(fPsi) IMPLICIT NONE - CLASS(meshCell3DCartTetra), INTENT(in):: self REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) @@ -338,127 +318,45 @@ MODULE moduleMesh3DCart END FUNCTION fPsiTetra !Derivative element function at coordinates Xi - PURE FUNCTION dPsiTetra(self, Xi, nNodes) RESULT(dPsi) + PURE FUNCTION dPsiTetra(Xi, nNodes) RESULT(dPsi) IMPLICIT NONE - CLASS(meshCell3DCartTetra), 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,:) = dPsiTetraXi1(Xi(2), Xi(3)) - dPsi(2,:) = dPsiTetraXi2(Xi(1), Xi(3)) - dPsi(3,:) = dPsiTetraXi3(Xi(1), Xi(2)) + dPsi(1,1:4) = (/ -1.D0, 1.D0, 0.D0, 0.D0 /) + dPsi(2,1:4) = (/ -1.D0, 0.D0, 1.D0, 0.D0 /) + dPsi(3,1:4) = (/ -1.D0, 0.D0, 0.D0, 1.D0 /) END FUNCTION dPsiTetra - !Derivative element function respect to Xi1 - PURE FUNCTION dPsiTetraXi1(Xi2, Xi3) RESULT(dPsiXi1) - IMPLICIT NONE - REAL(8), INTENT(in):: Xi2, Xi3 - REAL(8):: dPsiXi1(1:4) - - dPsiXi1(1) = -1.D0 - dPsiXi1(2) = 1.D0 - dPsiXi1(3) = 0.D0 - dPsiXi1(4) = 0.D0 - - END FUNCTION dPsiTetraXi1 - - !Derivative element function respect to Xi2 - PURE FUNCTION dPsiTetraXi2(Xi1, Xi3) RESULT(dPsiXi2) - IMPLICIT NONE - REAL(8), INTENT(in):: Xi1, Xi3 - REAL(8):: dPsiXi2(1:4) - - dPsiXi2(1) = -1.D0 - dPsiXi2(2) = 0.D0 - dPsiXi2(3) = 1.D0 - dPsiXi2(4) = 0.D0 - - END FUNCTION dPsiTetraXi2 - - !Derivative element function respect to Xi3 - PURE FUNCTION dPsiTetraXi3(Xi1, Xi2) RESULT(dPsiXi3) - IMPLICIT NONE - REAL(8), INTENT(in):: Xi1, Xi2 - REAL(8):: dPsiXi3(1:4) - - dPsiXi3(1) = -1.D0 - dPsiXi3(2) = 0.D0 - dPsiXi3(3) = 0.D0 - dPsiXi3(4) = 1.D0 - - END FUNCTION dPsiTetraXi3 - !Computes the derivatives in global coordinates - PURE SUBROUTINE partialDerTetra(self, nNodes, dPsi, dx, dy, dz) + PURE FUNCTION partialDerTetra(self, nNodes, dPsi) RESULT(pDer) IMPLICIT NONE CLASS(meshCell3DCartTetra), INTENT(in):: self INTEGER, INTENT(in):: nNodes REAL(8), INTENT(in):: dPsi(1:3, 1:nNodes) - REAL(8), INTENT(out), DIMENSION(1:3):: dx, dy, dz + REAL(8):: pDer(1:3, 1:3) - dx(1) = DOT_PRODUCT(dPsi(1,:), self%x) - dx(2) = DOT_PRODUCT(dPsi(2,:), self%x) - dx(3) = DOT_PRODUCT(dPsi(3,:), self%x) + pDer = 0.D0 - dy(1) = DOT_PRODUCT(dPsi(1,:), self%y) - dy(2) = DOT_PRODUCT(dPsi(2,:), self%y) - dy(3) = DOT_PRODUCT(dPsi(3,:), self%y) + pDer(1, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%x(1:4)), & + DOT_PRODUCT(dPsi(2,1:4), self%x(1:4)), & + DOT_PRODUCT(dPsi(3,1:4), self%x(1:4)) /) - dz(1) = DOT_PRODUCT(dPsi(1,:), self%z) - dz(2) = DOT_PRODUCT(dPsi(2,:), self%z) - dz(3) = DOT_PRODUCT(dPsi(3,:), self%z) + pDer(2, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%y(1:4)), & + DOT_PRODUCT(dPsi(2,1:4), self%y(1:4)), & + DOT_PRODUCT(dPsi(3,1:4), self%y(1:4)) /) - END SUBROUTINE partialDerTetra + pDer(3, 1:3) = (/ DOT_PRODUCT(dPsi(1,1:4), self%z(1:4)), & + DOT_PRODUCT(dPsi(2,1:4), self%z(1:4)), & + DOT_PRODUCT(dPsi(3,1:4), self%z(1:4)) /) - PURE FUNCTION elemKTetra(self, nNodes) RESULT(localK) - IMPLICIT NONE - - CLASS(meshCell3DCartTetra), 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 - - localK = 0.D0 - Xi = 0.D0 - !TODO: One point Gauss integral. Upgrade when possible - Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) - dPsi = self%dPsi(Xi, 4) - 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, nNodes, source) RESULT(localF) - IMPLICIT NONE - - CLASS(meshCell3DCartTetra), INTENT(in):: self - 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 - - localF = 0.D0 - Xi = 0.D0 - Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) - dPsi = self%dPsi(Xi, 4) - detJ = self%detJac(Xi, 4, dPsi) - fPsi = self%fPsi(Xi, 4) - f = DOT_PRODUCT(fPsi, source) - localF = f*fPsi*1.D0*detJ - - END FUNCTION elemFTetra + END FUNCTION partialDerTetra PURE FUNCTION gatherEFTetra(self, Xi) RESULT(array) IMPLICIT NONE @@ -502,6 +400,54 @@ MODULE moduleMesh3DCart END FUNCTION gatherMFTetra + PURE FUNCTION elemKTetra(self, nNodes) RESULT(localK) + IMPLICIT NONE + + CLASS(meshCell3DCartTetra), 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):: pDer(1:3, 1:3) + REAL(8):: invJ(1:3,1:3), detJ + + localK = 0.D0 + Xi = 0.D0 + !TODO: One point Gauss integral. Upgrade when possible + Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) + fPsi = self%fPsi(Xi, 4) + localK = MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*1.D0/detJ + + END FUNCTION elemKTetra + + PURE FUNCTION elemFTetra(self, nNodes, source) RESULT(localF) + IMPLICIT NONE + + CLASS(meshCell3DCartTetra), 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), dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) + REAL(8):: detJ, f + + localF = 0.D0 + Xi = 0.D0 + Xi = (/ 0.25D0, 0.25D0, 0.25D0 /) + 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 = f*fPsi*1.D0*detJ + + END FUNCTION elemFTetra + PURE FUNCTION insideTetra(Xi) RESULT(ins) IMPLICIT NONE @@ -515,121 +461,101 @@ MODULE moduleMesh3DCart END FUNCTION insideTetra - PURE FUNCTION getNodesTetra(self, nNodes) RESULT(n) - IMPLICIT NONE - - CLASS(meshCell3DCartTetra), 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 getNodesTetra - PURE FUNCTION phy2logTetra(self,r) RESULT(Xi) IMPLICIT NONE CLASS(meshCell3DCartTetra), INTENT(in):: self REAL(8), INTENT(in):: r(1:3) REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) REAL(8):: invJ(1:3, 1:3), detJ REAL(8):: deltaR(1:3) - REAL(8):: dPsi(1:3, 1:4) Xi = 0.D0 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, 4, dPsi) - detJ = self%detJac(Xi, 4, dPsi) + pDer = self%partialDer(4, dPsi) + invJ = self%invJac(pDer) + detJ = self%detJac(pDer) Xi = MATMUL(invJ, deltaR)/detJ END FUNCTION phy2logTetra - SUBROUTINE nextElementTetra(self, Xi, nextElement) + SUBROUTINE neighbourElementTetra(self, Xi, neighbourElement) IMPLICIT NONE CLASS(meshCell3DCartTetra), 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 !TODO: Review when connectivity XiArray = (/ Xi(3), 1.D0 - Xi(1) - Xi(2) - Xi(3), 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 CASE (4) - nextElement => self%e4 + neighbourElement => self%e4 END SELECT - END SUBROUTINE nextElementTetra + END SUBROUTINE neighbourElementTetra + + !Computes the element volume + PURE SUBROUTINE volumeTetra(self) + IMPLICIT NONE + + CLASS(meshCell3DCartTetra), INTENT(inout):: self + REAL(8):: Xi(1:3) + REAL(8):: dPsi(1:3, 1:4) + REAL(8):: pDer(1:3, 1:3) + + self%volume = 0.D0 + Xi = (/0.25D0, 0.25D0, 0.25D0/) + dPsi = self%dPsi(Xi, 4) + pDer = self%partialDer(4, dPsi) + self%volume = self%detJac(pDer) + + END SUBROUTINE volumeTetra !COMMON FUNCTIONS FOR CARTESIAN VOLUME ELEMENTS IN 3D !Computes element Jacobian determinant - PURE FUNCTION detJ3DCart(self, Xi, nNodes, dPsi_in) RESULT(dJ) + PURE FUNCTION detJ3DCart(pDer) RESULT(dJ) IMPLICIT NONE - CLASS(meshCell3DCart), 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:3), dy(1:3), dz(1:3) - IF (PRESENT(dPsi_in)) THEN - dPsi = dPsi_in - - ELSE - dPsi = self%dPsi(Xi, 4) - - END IF - - 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)) + dJ = pDer(1,1)*(pDer(2,2)*pDer(3,3) - pDer(2,3)*pDer(3,2)) & + - pDer(1,2)*(pDer(2,1)*pDer(3,3) - pDer(2,3)*pDer(3,1)) & + + pDer(1,3)*(pDer(2,1)*pDer(3,2) - pDer(2,2)*pDer(3,1)) END FUNCTION detJ3DCart - PURE FUNCTION invJ3DCart(self, Xi, nNodes, dPsi_in) RESULT(invJ) + PURE FUNCTION invJ3DCart(pDer) RESULT(invJ) IMPLICIT NONE - CLASS(meshCell3DCart), 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), DIMENSION(1:3):: dx, dy, dz + REAL(8), INTENT(in):: pDer(1:3, 1:3) REAL(8):: invJ(1:3,1:3) - IF(PRESENT(dPsi_in)) THEN - dPsi=dPsi_in + invJ(1,1:3) = (/ (pDer(2,2)*pDer(3,3) - pDer(2,3)*pDer(3,2)), & + -(pDer(2,1)*pDer(3,3) - pDer(2,3)*pDer(3,1)), & + (pDer(2,1)*pDer(3,2) - pDer(2,2)*pDer(3,1)) /) - ELSE - dPsi = self%dPsi(Xi, 4) + invJ(2,1:3) = (/ -(pDer(1,2)*pDer(3,3) - pDer(1,3)*pDer(3,2)), & + (pDer(1,1)*pDer(3,3) - pDer(1,3)*pDer(3,1)), & + -(pDer(1,1)*pDer(3,2) - pDer(1,2)*pDer(3,1)) /) - END IF - - 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)) - - invJ(2,1) = -(dx(2)*dz(3) - dx(3)*dz(2)) - invJ(2,2) = (dx(1)*dz(3) - dx(3)*dz(1)) - invJ(2,3) = -(dx(1)*dz(2) - dx(2)*dz(1)) - - invJ(3,1) = (dx(2)*dy(3) - dx(3)*dy(2)) - invJ(3,2) = -(dx(1)*dy(3) - dx(3)*dy(1)) - invJ(3,3) = (dx(1)*dy(2) - dx(2)*dy(1)) + invJ(3,1:3) = (/ (pDer(1,2)*pDer(2,3) - pDer(1,3)*pDer(2,2)), & + -(pDer(1,1)*pDer(2,3) - pDer(1,3)*pDer(2,1)), & + (pDer(1,1)*pDer(2,2) - pDer(1,2)*pDer(2,1)) /) invJ = TRANSPOSE(invJ) diff --git a/src/modules/mesh/moduleMesh.f90 b/src/modules/mesh/moduleMesh.f90 index 6561850..97ce691 100644 --- a/src/modules/mesh/moduleMesh.f90 +++ b/src/modules/mesh/moduleMesh.f90 @@ -24,8 +24,10 @@ MODULE moduleMesh !Lock indicator for scattering INTEGER(KIND=OMP_LOCK_KIND):: lock CONTAINS + !DEFERED PROCEDURES PROCEDURE(initNode_interface), DEFERRED, PASS:: init PROCEDURE(getCoord_interface), DEFERRED, PASS:: getCoordinates + !GENERIC PROCEDURES PROCEDURE, PASS:: resetOutput END TYPE meshNode @@ -83,6 +85,7 @@ MODULE moduleMesh !Physical surface for the edge INTEGER:: physicalSurface CONTAINS + !DEFERED PROCEDURES PROCEDURE(initEdge_interface), DEFERRED, PASS:: init PROCEDURE(getNodesEdge_interface), DEFERRED, PASS:: getNodes PROCEDURE(intersectionEdge_interface), DEFERRED, PASS:: intersection @@ -166,37 +169,41 @@ MODULE moduleMesh !Total weight of particles inside cell REAL(8), ALLOCATABLE:: totalWeight(:) CONTAINS + !DEFERRED PROCEDURES !Init the cell - PROCEDURE(initCell_interface), DEFERRED, PASS:: init + PROCEDURE(initCell_interface), DEFERRED, PASS:: init !Get the index of the nodes - PROCEDURE(getNodesCell_interface), DEFERRED, PASS:: getNodes + PROCEDURE(getNodesCell_interface), DEFERRED, PASS:: getNodes !Calculate random position on the cell - PROCEDURE(randPosVol_interface), DEFERRED, PASS:: randPos + PROCEDURE(randPosCell_interface), DEFERRED, PASS:: randPos !Obtain functions and values of cell natural functions - PROCEDURE(fPsi_interface), DEFERRED, PASS:: fPsi - PROCEDURE(dPsi_interface), DEFERRED, PASS:: dPsi - PROCEDURE(detJac_interface), DEFERRED, PASS:: detJac - PROCEDURE(invJac_interface), DEFERRED, PASS:: invJac + PROCEDURE(fPsi_interface), DEFERRED, NOPASS:: fPsi + PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi + PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer + PROCEDURE(detJac_interface), DEFERRED, NOPASS:: detJac + PROCEDURE(invJac_interface), DEFERRED, NOPASS:: invJac + !Procedures to get specific values in the node + PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherElectricField + PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherMagneticField + !Compute K and F to solve PDE on the mesh + PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK + PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF + !Check if particle is inside the cell + PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside + !Convert physical coordinates (r) into logical coordinates (Xi) + PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log + !Returns the neighbour element based on particle position outside the cell + PROCEDURE(neighbourElement_interface), DEFERRED, PASS:: neighbourElement !Scatter properties of particles on cell nodes PROCEDURE, PASS:: scatter + !Subroutine to find in which cell a particle is located + PROCEDURE, PASS:: findCell !Gather value and spatial derivative on the nodes at position Xi PROCEDURE, PASS, PRIVATE:: gatherF_scalar PROCEDURE, PASS, PRIVATE:: gatherF_array PROCEDURE, PASS, PRIVATE:: gatherDF_scalar GENERIC:: gatherF => gatherF_scalar, gatherF_array GENERIC:: gatherDF => gatherDF_scalar - !Procedures to get specific values in the node - PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherElectricField - PROCEDURE(gatherArray_interface), DEFERRED, PASS:: gatherMagneticField - !Compute K and F to solve PDE on the mesh - PROCEDURE(elemK_interface), DEFERRED, PASS:: elemK - PROCEDURE(elemF_interface), DEFERRED, PASS:: elemF - !Subroutines to find in which cell a particle is located - PROCEDURE, PASS:: findCell - PROCEDURE(inside_interface), DEFERRED, NOPASS:: inside - PROCEDURE(nextElement_interface), DEFERRED, PASS:: nextElement - !Convert physical coordinates (r) into logical coordinates (Xi) - PROCEDURE(phy2log_interface), DEFERRED, PASS:: phy2log END TYPE meshCell @@ -219,40 +226,44 @@ MODULE moduleMesh END FUNCTION getNodesCell_interface - PURE FUNCTION fPsi_interface(self, Xi, nNodes) RESULT(fPsi) + FUNCTION randPosCell_interface(self) RESULT(r) IMPORT:: meshCell CLASS(meshCell), INTENT(in):: self + REAL(8):: r(1:3) + + END FUNCTION randPosCell_interface + + PURE FUNCTION fPsi_interface(Xi, nNodes) RESULT(fPsi) REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: fPsi(1:nNodes) END FUNCTION fPsi_interface - PURE FUNCTION dPsi_interface(self, Xi, nNodes) RESULT(dPsi) - IMPORT:: meshCell - CLASS(meshCell), INTENT(in):: self + PURE FUNCTION dPsi_interface(Xi, nNodes) RESULT(dPsi) REAL(8), INTENT(in):: Xi(1:3) INTEGER, INTENT(in):: nNodes REAL(8):: dPsi(1:3, 1:nNodes) END FUNCTION dPsi_interface - PURE FUNCTION detJac_interface(self, Xi, nNodes, dPsi_in) RESULT(dJ) + PURE FUNCTION partialDer_interface(self, nNodes, dPsi) RESULT(pDer) IMPORT:: meshCell CLASS(meshCell), 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):: dPsi(1:3, 1:nNodes) + REAL(8):: pDer(1:3, 1:3) + + END FUNCTION partialDer_interface + + PURE FUNCTION detJac_interface(pDer) RESULT(dJ) + REAL(8), INTENT(in):: pDer(1:3,1:3) REAL(8):: dJ END FUNCTION detJac_interface - 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) - INTEGER, INTENT(in):: nNodes - REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:3,1:nNodes) + PURE FUNCTION invJac_interface(pDer) RESULT(invJ) + REAL(8), INTENT(in):: pDer(1:3,1:3) REAL(8):: invJ(1:3,1:3) END FUNCTION invJac_interface @@ -282,13 +293,12 @@ MODULE moduleMesh END FUNCTION elemF_interface - SUBROUTINE nextElement_interface(self, Xi, nextElement) - IMPORT:: meshCell, meshElement - CLASS(meshCell), INTENT(in):: self + PURE FUNCTION inside_interface(Xi) RESULT(ins) + IMPORT:: meshCell REAL(8), INTENT(in):: Xi(1:3) - CLASS(meshElement), POINTER, INTENT(out):: nextElement + LOGICAL:: ins - END SUBROUTINE nextElement_interface + END FUNCTION inside_interface PURE FUNCTION phy2log_interface(self,r) RESULT(Xi) IMPORT:: meshCell @@ -298,19 +308,13 @@ MODULE moduleMesh END FUNCTION phy2log_interface - PURE FUNCTION inside_interface(Xi) RESULT(ins) - IMPORT:: meshCell - REAL(8), INTENT(in):: Xi(1:3) - LOGICAL:: ins - - END FUNCTION inside_interface - - FUNCTION randPosVol_interface(self) RESULT(r) - IMPORT:: meshCell + SUBROUTINE neighbourElement_interface(self, Xi, neighbourElement) + IMPORT:: meshCell, meshElement CLASS(meshCell), INTENT(in):: self - REAL(8):: r(1:3) + REAL(8), INTENT(in):: Xi(1:3) + CLASS(meshElement), POINTER, INTENT(out):: neighbourElement - END FUNCTION randPosVol_interface + END SUBROUTINE neighbourElement_interface END INTERFACE @@ -332,11 +336,13 @@ MODULE moduleMesh TYPE(meshNodeCont), ALLOCATABLE:: nodes(:) !Array of cell elements TYPE(meshCellCont), ALLOCATABLE:: cells(:) + !PROCEDURES SPECIFIC OF FILE TYPE PROCEDURE(readMesh_interface), POINTER, PASS:: readMesh => NULL() PROCEDURE(readInitial_interface), POINTER, NOPASS:: readInitial => NULL() PROCEDURE(connectMesh_interface), POINTER, PASS:: connectMesh => NULL() PROCEDURE(printColl_interface), POINTER, PASS:: printColl => NULL() CONTAINS + !GENERIC PROCEDURES PROCEDURE, PASS:: doCollisions END TYPE meshGeneric @@ -345,7 +351,6 @@ MODULE moduleMesh !Reads the mesh from a file SUBROUTINE readMesh_interface(self, filename) IMPORT meshGeneric - CLASS(meshGeneric), INTENT(inout):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: filename @@ -363,7 +368,6 @@ MODULE moduleMesh !Connects cell and edges to the mesh SUBROUTINE connectMesh_interface(self) IMPORT meshGeneric - CLASS(meshGeneric), INTENT(inout):: self END SUBROUTINE connectMesh_interface @@ -371,7 +375,6 @@ MODULE moduleMesh !Prints number of collisions in each cell SUBROUTINE printColl_interface(self, t) IMPORT meshGeneric - CLASS(meshGeneric), INTENT(inout):: self INTEGER, INTENT(in):: t @@ -388,28 +391,21 @@ MODULE moduleMesh REAL(8), ALLOCATABLE, DIMENSION(:,:):: K !Permutation matrix for P L U factorization INTEGER, ALLOCATABLE, DIMENSION(:,:):: IPIV + !PROCEDURES SPECIFIC OF FILE TYPE PROCEDURE(printOutput_interface), POINTER, PASS:: printOutput => NULL() PROCEDURE(printEM_interface), POINTER, PASS:: printEM => NULL() PROCEDURE(doCoulomb_interface), POINTER, PASS:: doCoulomb => NULL() PROCEDURE(printAverage_interface), POINTER, PASS:: printAverage => NULL() CONTAINS + !GENERIC PROCEDURES PROCEDURE, PASS:: constructGlobalK END TYPE meshParticles ABSTRACT INTERFACE - !Perform Coulomb Scattering - SUBROUTINE doCoulomb_interface(self) - IMPORT meshParticles - - CLASS(meshParticles), INTENT(inout):: self - - END SUBROUTINE doCoulomb_interface - !Prints Species data SUBROUTINE printOutput_interface(self, t) IMPORT meshParticles - CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: t @@ -418,21 +414,25 @@ MODULE moduleMesh !Prints EM info SUBROUTINE printEM_interface(self, t) IMPORT meshParticles - CLASS(meshParticles), INTENT(in):: self INTEGER, INTENT(in):: t END SUBROUTINE printEM_interface + !Perform Coulomb Scattering + SUBROUTINE doCoulomb_interface(self) + IMPORT meshParticles + CLASS(meshParticles), INTENT(inout):: self + + END SUBROUTINE doCoulomb_interface + !Prints average values SUBROUTINE printAverage_interface(self) IMPORT meshParticles - CLASS(meshParticles), INTENT(in):: self END SUBROUTINE printAverage_interface - END INTERFACE TYPE(meshParticles), TARGET:: mesh @@ -440,6 +440,7 @@ MODULE moduleMesh !Collision (MCC) mesh TYPE, EXTENDS(meshGeneric):: meshCollisions CONTAINS + !GENERIC PROCEDURES END TYPE meshCollisions @@ -448,7 +449,6 @@ MODULE moduleMesh ABSTRACT INTERFACE SUBROUTINE readMeshColl_interface(self, filename) IMPORT meshCollisions - CLASS(meshCollisions), INTENT(inout):: self CHARACTER(:), ALLOCATABLE, INTENT(in):: filename @@ -577,12 +577,14 @@ MODULE moduleMesh REAL(8), INTENT(in):: valNodes(1:nNodes) REAL(8):: df(1:3) REAL(8):: dPsi(1:3, 1:nNodes) + REAL(8):: pDer(1:3,1:3) REAL(8):: dPsiR(1:3, 1:nNodes) REAL(8):: invJ(1:3, 1:3), detJ dPsi = self%dPsi(Xi, nNodes) - detJ = self%detJac(Xi, nNodes, dPsi) - invJ = self%invJac(Xi, nNodes, dPsi) + pDer = self%partialDer(nNodes, dPsi) + detJ = self%detJac(pDer) + invJ = self%invJac(pDer) dPsiR = MATMUL(invJ, dPsi)/detJ df = (/ DOT_PRODUCT(dPsiR(1,:), valNodes), & DOT_PRODUCT(dPsiR(2,:), valNodes), & @@ -637,7 +639,7 @@ MODULE moduleMesh CLASS(particle), INTENT(inout), TARGET:: part CLASS(meshCell), OPTIONAL, INTENT(in):: oldCell REAL(8):: Xi(1:3) - CLASS(meshElement), POINTER:: nextElement + CLASS(meshElement), POINTER:: neighbourElement INTEGER:: sp Xi = self%phy2log(part%r) @@ -655,16 +657,16 @@ MODULE moduleMesh ELSE !If not, searches for a neighbour and repeats the process. - CALL self%nextElement(Xi, nextElement) + CALL self%neighbourElement(Xi, neighbourElement) !Defines the next step - SELECT TYPE(nextElement) + SELECT TYPE(neighbourElement) CLASS IS(meshCell) !Particle moved to new cell, repeat find procedure - CALL nextElement%findCell(part, self) + CALL neighbourElement%findCell(part, self) CLASS IS (meshEdge) !Particle encountered a surface, apply boundary - CALL nextElement%fBoundary(part%species%n)%apply(nextElement,part) + CALL neighbourElement%fBoundary(part%species%n)%apply(neighbourElement,part) !If particle is still inside the domain, call findCell IF (part%n_in) THEN @@ -709,7 +711,7 @@ MODULE moduleMesh LOGICAL:: found CLASS(meshCell), POINTER:: cell REAL(8), DIMENSION(1:3):: Xi - CLASS(meshElement), POINTER:: nextElement + CLASS(meshElement), POINTER:: neighbourElement INTEGER:: sp found = .FALSE. @@ -727,11 +729,11 @@ MODULE moduleMesh found = .TRUE. ELSE - CALL cell%nextElement(Xi, nextElement) - SELECT TYPE(nextElement) + CALL cell%neighbourElement(Xi, neighbourElement) + SELECT TYPE(neighbourElement) CLASS IS(meshCell) !Try next element - cell => nextElement + cell => neighbourElement CLASS DEFAULT !Should never happend, but just in case, stops loops