JGonzalez
386ddd82dd
Probes are now written at the 0 iteration. Additionally, and this shouldn't be done, some small changes to the quad elements. This should be done in a separate commit, but I'm lazy.
1471 lines
40 KiB
Fortran
1471 lines
40 KiB
Fortran
!moduleMesh2DCyl: 2D aXial symmetric extension of generic mesh from GMSH format.
|
|
! x == z
|
|
! y == r
|
|
! z == theta (unused)
|
|
MODULE moduleMesh2DCyl
|
|
USE moduleMesh
|
|
USE moduleMeshBoundary
|
|
IMPLICIT NONE
|
|
|
|
!Values for Gauss integral
|
|
REAL(8), PARAMETER:: corQuad(1:3) = (/ -DSQRT(3.D0/5.D0), 0.D0, DSQRT(3.D0/5.D0) /)
|
|
REAL(8), PARAMETER:: wQuad(1:3) = (/ 5.D0/9.D0, 8.D0/9.D0, 5.D0/9.D0 /)
|
|
|
|
REAL(8), PARAMETER:: Xi1Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 3.D0/5.D0, 1.D0/5.D0 /)
|
|
REAL(8), PARAMETER:: Xi2Tria(1:4) = (/ 1.D0/3.D0, 1.D0/5.D0, 1.D0/5.D0, 3.D0/5.D0 /)
|
|
REAL(8), PARAMETER:: wTria(1:4) = (/ -27.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0, 25.D0/96.D0 /)
|
|
|
|
TYPE, PUBLIC, EXTENDS(meshNode):: meshNode2DCyl
|
|
!Element coordinates
|
|
REAL(8):: r = 0.D0, z = 0.D0
|
|
CONTAINS
|
|
!meshNode DEFERRED PROCEDURES
|
|
PROCEDURE, PASS:: init => initNode2DCyl
|
|
PROCEDURE, PASS:: getCoordinates => getCoord2DCyl
|
|
|
|
END TYPE meshNode2DCyl
|
|
|
|
TYPE, PUBLIC, EXTENDS(meshEdge):: meshEdge2DCyl
|
|
!Element coordinates
|
|
REAL(8):: r(1:2) = 0.D0, z(1:2) = 0.D0
|
|
!Connectivity to nodes
|
|
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL()
|
|
CONTAINS
|
|
!meshEdge DEFERRED PROCEDURES
|
|
PROCEDURE, PASS:: init => initEdge2DCyl
|
|
PROCEDURE, PASS:: getNodes => getNodes2DCyl
|
|
PROCEDURE, PASS:: intersection => intersection2DCylEdge
|
|
PROCEDURE, PASS:: randPos => randPosEdge
|
|
|
|
END TYPE meshEdge2DCyl
|
|
|
|
!Quadrilateral volume element
|
|
TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylQuad
|
|
!Element coordinates
|
|
REAL(8):: r(1:4) = 0.D0, z(1:4) = 0.D0
|
|
!Connectivity to nodes
|
|
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL(), n4 => NULL()
|
|
!Connectivity to adjacent elements
|
|
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL(), e4 => NULL()
|
|
|
|
CONTAINS
|
|
!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, PRIVATE:: calculateVolume => volumeQuad
|
|
|
|
END TYPE meshCell2DCylQuad
|
|
|
|
!Triangular volume element
|
|
TYPE, PUBLIC, EXTENDS(meshCell):: meshCell2DCylTria
|
|
!Element coordinates
|
|
REAL(8):: r(1:3) = 0.D0, z(1:3) = 0.D0
|
|
!Connectivity to nodes
|
|
CLASS(meshNode), POINTER:: n1 => NULL(), n2 => NULL(), n3 => NULL()
|
|
!Connectivity to adjacent elements
|
|
CLASS(meshElement), POINTER:: e1 => NULL(), e2 => NULL(), e3 => NULL()
|
|
|
|
CONTAINS
|
|
!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, PRIVATE:: calculateVolume => volumeTria
|
|
|
|
END TYPE meshCell2DCylTria
|
|
|
|
CONTAINS
|
|
!NODE FUNCTIONS
|
|
!Init node element
|
|
SUBROUTINE initNode2DCyl(self, n, r)
|
|
USE moduleSpecies
|
|
USE moduleRefParam
|
|
USE OMP_LIB
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshNode2DCyl), INTENT(out):: self
|
|
INTEGER, INTENT(in):: n
|
|
REAL(8), INTENT(in):: r(1:3)
|
|
|
|
self%n = n
|
|
self%z = r(1)/L_ref
|
|
self%r = r(2)/L_ref
|
|
!Node volume, to be determined in mesh
|
|
self%v = 0.D0
|
|
|
|
!Allocates output:
|
|
ALLOCATE(self%output(1:nSpecies))
|
|
|
|
CALL OMP_INIT_LOCK(self%lock)
|
|
|
|
END SUBROUTINE initNode2DCyl
|
|
|
|
!Get coordinates from node
|
|
PURE FUNCTION getCoord2DCyl(self) RESULT(r)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshNode2DCyl), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
|
|
r = (/self%z, self%r, 0.D0/)
|
|
|
|
END FUNCTION getCoord2DCyl
|
|
|
|
!EDGE FUNCTIONS
|
|
!Init edge element
|
|
SUBROUTINE initEdge2DCyl(self, n, p, bt, physicalSurface)
|
|
USE moduleSpecies
|
|
USE moduleBoundary
|
|
USE moduleErrors
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshEdge2DCyl), INTENT(out):: self
|
|
INTEGER, INTENT(in):: n
|
|
INTEGER, INTENT(in):: p(:)
|
|
INTEGER, INTENT(in):: bt
|
|
INTEGER, INTENT(in):: physicalSurface
|
|
REAL(8), DIMENSION(1:3):: r1, r2
|
|
INTEGER:: s
|
|
|
|
self%n = n
|
|
self%nNodes = SIZE(p)
|
|
self%n1 => mesh%nodes(p(1))%obj
|
|
self%n2 => mesh%nodes(p(2))%obj
|
|
!Get element coordinates
|
|
r1 = self%n1%getCoordinates()
|
|
r2 = self%n2%getCoordinates()
|
|
self%z = (/r1(1), r2(1)/)
|
|
self%r = (/r1(2), r2(2)/)
|
|
self%weight = SUM(self%r)*5.D-1
|
|
!Normal vector
|
|
self%normal = (/ -(self%r(2)-self%r(1)), &
|
|
self%z(2)-self%z(1) , &
|
|
0.D0 /)
|
|
self%normal = self%normal/NORM2(self%normal)
|
|
|
|
!Boundary index
|
|
self%boundary => boundary(bt)
|
|
ALLOCATE(self%fboundary(1:nSpecies))
|
|
!Assign functions to boundary
|
|
DO s = 1, nSpecies
|
|
CALL pointBoundaryFunction(self, s)
|
|
|
|
END DO
|
|
|
|
!Physical surface
|
|
self%physicalSurface = physicalSurface
|
|
|
|
END SUBROUTINE initEdge2DCyl
|
|
|
|
!Get nodes from edge
|
|
PURE FUNCTION getNodes2DCyl(self, nNodes) RESULT(n)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshEdge2DCyl), INTENT(in):: self
|
|
INTEGER, INTENT(in):: nNodes
|
|
INTEGER:: n(1:nNodes)
|
|
|
|
n = (/self%n1%n, self%n2%n /)
|
|
|
|
END FUNCTION getNodes2DCyl
|
|
|
|
!Calculate intersection between position and edge
|
|
PURE FUNCTION intersection2DCylEdge(self, r0) RESULT(r)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshEdge2DCyl), INTENT(in):: self
|
|
REAL(8), DIMENSION(1:3), INTENT(in):: r0
|
|
REAL(8), DIMENSION(1:3):: r
|
|
REAL(8), DIMENSION(1:3):: edge0, edgeV
|
|
REAL(8):: tI
|
|
|
|
edge0 = (/self%z(1), self%r(1), 0.D0 /)
|
|
edgeV = (/self%z(2), self%r(2), 0.D0 /) - edge0
|
|
|
|
tI = DOT_PRODUCT(r0 - edge0, edgeV)/DOT_PRODUCT(edgeV, edgeV)
|
|
|
|
r = edge0 + tI*edgeV
|
|
|
|
END FUNCTION intersection2DCylEdge
|
|
|
|
!Calculate a random position in edge
|
|
FUNCTION randPosEdge(self) RESULT(r)
|
|
USE moduleRandom
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshEdge2DCyl), INTENT(in):: self
|
|
REAL(8):: rnd
|
|
REAL(8):: r(1:3)
|
|
REAL(8):: dr, dz
|
|
|
|
rnd = random()
|
|
dr = self%r(2) - self%r(1)
|
|
dz = self%z(2) - self%z(1)
|
|
IF (dr /= 0.D0) THEN
|
|
r(2) = dr * DSQRT(rnd) + self%r(1)
|
|
r(1) = dz * (r(2) - self%r(1))/dr + self%z(1)
|
|
|
|
ELSE
|
|
r(2) = self%r(1)
|
|
r(1) = dz * rnd + self%z(1)
|
|
|
|
END IF
|
|
|
|
r(3) = 0.D0
|
|
|
|
END FUNCTION randPosEdge
|
|
|
|
!VOLUME FUNCTIONS
|
|
!QUAD FUNCTIONS
|
|
!Init element
|
|
SUBROUTINE initCellQuad2DCyl(self, n, p, nodes)
|
|
USE moduleRefParam
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylQuad), INTENT(out):: self
|
|
INTEGER, INTENT(in):: n
|
|
INTEGER, INTENT(in):: p(:)
|
|
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
|
|
self%n4 => nodes(p(4))%obj
|
|
!Get element coordinates
|
|
r1 = self%n1%getCoordinates()
|
|
r2 = self%n2%getCoordinates()
|
|
r3 = self%n3%getCoordinates()
|
|
r4 = self%n4%getCoordinates()
|
|
self%z = (/r1(1), r2(1), r3(1), r4(1)/)
|
|
self%r = (/r1(2), r2(2), r3(2), r4(2)/)
|
|
|
|
!Assign node volume
|
|
CALL self%calculateVolume()
|
|
|
|
CALL OMP_INIT_LOCK(self%lock)
|
|
|
|
ALLOCATE(self%listPart_in(1:nSpecies))
|
|
ALLOCATE(self%totalWeight(1:nSpecies))
|
|
|
|
END SUBROUTINE initCellQuad2DCyl
|
|
|
|
!Get 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
|
|
|
|
!Random position in quadrilateral volume
|
|
FUNCTION randPosCellQuad(self) RESULT(r)
|
|
USE moduleRandom
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylQuad), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
REAL(8):: Xi(1:3)
|
|
REAL(8):: fPsi(1:4)
|
|
|
|
Xi(1) = random(-1.D0, 1.D0)
|
|
Xi(2) = random(-1.D0, 1.D0)
|
|
Xi(3) = 0.D0
|
|
|
|
fPsi = self%fPsi(Xi, 4)
|
|
|
|
r(1) = DOT_PRODUCT(fPsi, self%z)
|
|
r(2) = DOT_PRODUCT(fPsi, self%r)
|
|
r(3) = 0.D0
|
|
|
|
END FUNCTION randPosCellQuad
|
|
|
|
!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(meshCell2DCylQuad), INTENT(in):: self
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8), INTENT(in):: dPsi(1:3,1:nNodes)
|
|
REAL(8):: pDer(1:3, 1:3)
|
|
|
|
pDer = 0.D0
|
|
|
|
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)) /)
|
|
pDer(3,3) = 1.D0
|
|
|
|
END FUNCTION partialDerQuad
|
|
|
|
PURE FUNCTION gatherEFQuad(self, Xi) RESULT(array)
|
|
IMPLICIT NONE
|
|
CLASS(meshCell2DCylQuad), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
REAL(8):: array(1:3)
|
|
REAL(8):: phi(1:4)
|
|
|
|
phi = (/ self%n1%emData%phi, &
|
|
self%n2%emData%phi, &
|
|
self%n3%emData%phi, &
|
|
self%n4%emData%phi /)
|
|
|
|
array = -self%gatherDF(Xi, 4, phi)
|
|
|
|
END FUNCTION gatherEFQuad
|
|
|
|
PURE FUNCTION gatherMFQuad(self, Xi) RESULT(array)
|
|
IMPLICIT NONE
|
|
CLASS(meshCell2DCylQuad), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
REAL(8):: array(1:3)
|
|
REAL(8):: B(1:4,1:3)
|
|
|
|
B(:,1) = (/ self%n1%emData%B(1), &
|
|
self%n2%emData%B(1), &
|
|
self%n3%emData%B(1), &
|
|
self%n4%emData%B(1) /)
|
|
|
|
B(:,2) = (/ self%n1%emData%B(2), &
|
|
self%n2%emData%B(2), &
|
|
self%n3%emData%B(2), &
|
|
self%n4%emData%B(2) /)
|
|
|
|
B(:,3) = (/ self%n1%emData%B(3), &
|
|
self%n2%emData%B(3), &
|
|
self%n3%emData%B(3), &
|
|
self%n4%emData%B(3) /)
|
|
|
|
array = self%gatherF(Xi, 4, B)
|
|
|
|
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):: 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)
|
|
r = self%gatherF(Xi, 4, 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)
|
|
REAL(8):: 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 Xi is inside the element
|
|
PURE FUNCTION insideQuad(Xi) RESULT(ins)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
LOGICAL:: ins
|
|
|
|
ins = (Xi(1) >= -1.D0 .AND. Xi(1) <= 1.D0) .AND. &
|
|
(Xi(2) >= -1.D0 .AND. Xi(2) <= 1.D0)
|
|
|
|
END FUNCTION insideQuad
|
|
|
|
!Transform physical coordinates to element coordinates
|
|
PURE FUNCTION phy2logQuad(self,r) RESULT(Xi)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylQuad), INTENT(in):: self
|
|
REAL(8), INTENT(in):: r(1:3)
|
|
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
|
|
conv = 1.D0
|
|
XiO = 0.D0
|
|
|
|
DO WHILE(conv > 1.D-4)
|
|
dPsi = self%dPsi(XiO, 4)
|
|
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), &
|
|
0.D0 /) - r
|
|
Xi = XiO - MATMUL(invJ, f)/detJ
|
|
conv = MAXVAL(DABS(Xi-XiO),1)
|
|
XiO = Xi
|
|
|
|
END DO
|
|
|
|
END FUNCTION phy2logQuad
|
|
|
|
!Get the neighbour element for a logical position Xi
|
|
SUBROUTINE neighbourElementQuad(self, Xi, neighbourElement)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylQuad), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
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(neighbourElement)
|
|
SELECT CASE (nextInt)
|
|
CASE (1)
|
|
neighbourElement => self%e1
|
|
CASE (2)
|
|
neighbourElement => self%e2
|
|
CASE (3)
|
|
neighbourElement => self%e3
|
|
CASE (4)
|
|
neighbourElement => self%e4
|
|
END SELECT
|
|
|
|
END SUBROUTINE neighbourElementQuad
|
|
|
|
!Compute element volume
|
|
PURE SUBROUTINE volumeQuad(self)
|
|
USE moduleConstParam, ONLY: PI8
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylQuad), INTENT(inout):: self
|
|
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)
|
|
|
|
self%volume = 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)
|
|
r = DOT_PRODUCT(fPsi,self%r)
|
|
!Computes total volume of the cell
|
|
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%n1%v = self%n1%v + fPsi(1)*r*detJ*PI8
|
|
Xi = (/ 5.D-1, -5.D-1, 0.D0/)
|
|
r = self%gatherF(Xi, 4, self%r)
|
|
self%n2%v = self%n2%v + fPsi(2)*r*detJ*PI8
|
|
Xi = (/ 5.D-1, 5.D-1, 0.D0/)
|
|
r = self%gatherF(Xi, 4, self%r)
|
|
self%n3%v = self%n3%v + fPsi(3)*r*detJ*PI8
|
|
Xi = (/-5.D-1, 5.D-1, 0.D0/)
|
|
r = self%gatherF(Xi, 4, self%r)
|
|
self%n4%v = self%n4%v + fPsi(4)*r*detJ*PI8
|
|
|
|
END SUBROUTINE volumeQuad
|
|
|
|
!TRIA FUNCTIONS
|
|
!Init element
|
|
SUBROUTINE initCellTria2DCyl(self, n, p, nodes)
|
|
USE moduleRefParam
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylTria), INTENT(out):: self
|
|
INTEGER, INTENT(in):: n
|
|
INTEGER, INTENT(in):: p(:)
|
|
TYPE(meshNodeCont), INTENT(in), TARGET:: nodes(:)
|
|
REAL(8), DIMENSION(1:3):: r1, r2, r3
|
|
|
|
!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
|
|
!Get element coordinates
|
|
r1 = self%n1%getCoordinates()
|
|
r2 = self%n2%getCoordinates()
|
|
r3 = self%n3%getCoordinates()
|
|
self%z = (/r1(1), r2(1), r3(1)/)
|
|
self%r = (/r1(2), r2(2), r3(2)/)
|
|
!Assign node volume
|
|
CALL self%calculateVolume()
|
|
|
|
CALL OMP_INIT_LOCK(self%lock)
|
|
|
|
ALLOCATE(self%listPart_in(1:nSpecies))
|
|
ALLOCATE(self%totalWeight(1:nSpecies))
|
|
|
|
END SUBROUTINE initCellTria2DCyl
|
|
|
|
!Random position in cell
|
|
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 cell
|
|
FUNCTION randPosCellTria(self) RESULT(r)
|
|
USE moduleRandom
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylTria), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
REAL(8):: Xi(1:3)
|
|
REAL(8):: fPsi(1:3)
|
|
|
|
Xi(1) = random( 0.D0, 1.D0)
|
|
Xi(2) = random( 0.D0, 1.D0 - Xi(1))
|
|
Xi(3) = 0.D0
|
|
|
|
fPsi = self%fPsi(Xi, 4)
|
|
|
|
r(1) = DOT_PRODUCT(fPsi, self%z)
|
|
r(2) = DOT_PRODUCT(fPsi, self%r)
|
|
r(3) = 0.D0
|
|
|
|
END FUNCTION randPosCellTria
|
|
|
|
!Compute element functions in point Xi
|
|
PURE FUNCTION fPsiTria(Xi, nNodes) RESULT(fPsi)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
INTEGER, INTENT(in):: nNodes
|
|
REAL(8):: fPsi(1:nNodes)
|
|
|
|
fPsi(1) = 1.D0 - Xi(1) - Xi(2)
|
|
fPsi(2) = Xi(1)
|
|
fPsi(3) = Xi(2)
|
|
|
|
END FUNCTION fPsiTria
|
|
|
|
!Compute element derivative functions in point Xi
|
|
PURE FUNCTION dPsiTria(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, 1.D0, 0.D0 /)
|
|
dPsi(2,:) = (/ -1.D0, 0.D0, 1.D0 /)
|
|
|
|
END FUNCTION dPsiTria
|
|
|
|
!Compute the derivatives in global coordinates
|
|
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):: pDer(1:3, 1:3)
|
|
|
|
pDer = 0.D0
|
|
|
|
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)) /)
|
|
|
|
END FUNCTION partialDerTria
|
|
|
|
!Gather electric field at position Xi
|
|
PURE FUNCTION gatherEFTria(self, Xi) RESULT(array)
|
|
IMPLICIT NONE
|
|
CLASS(meshCell2DCylTria), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
REAL(8):: array(1:3)
|
|
REAL(8):: phi(1:3)
|
|
|
|
phi = (/ self%n1%emData%phi, &
|
|
self%n2%emData%phi, &
|
|
self%n3%emData%phi /)
|
|
|
|
array = -self%gatherDF(Xi, 3, phi)
|
|
|
|
END FUNCTION gatherEFTria
|
|
|
|
!Gather magnetic field at position Xi
|
|
PURE FUNCTION gatherMFTria(self, Xi) RESULT(array)
|
|
IMPLICIT NONE
|
|
CLASS(meshCell2DCylTria), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
REAL(8):: array(1:3)
|
|
REAL(8):: B(1:3,1:3)
|
|
|
|
B(:,1) = (/ self%n1%emData%B(1), &
|
|
self%n2%emData%B(1), &
|
|
self%n3%emData%B(1) /)
|
|
|
|
B(:,2) = (/ self%n1%emData%B(2), &
|
|
self%n2%emData%B(2), &
|
|
self%n3%emData%B(2) /)
|
|
|
|
B(:,3) = (/ self%n1%emData%B(3), &
|
|
self%n2%emData%B(3), &
|
|
self%n3%emData%B(3) /)
|
|
|
|
array = self%gatherF(Xi, 3, B)
|
|
|
|
END FUNCTION gatherMFTria
|
|
|
|
!Compute cell 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):: dPsi(1:3,1:3)
|
|
REAL(8):: r
|
|
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)
|
|
r = self%gatherF(Xi, 3, self%r)
|
|
localK = localK + MATMUL(TRANSPOSE(MATMUL(invJ,dPsi)),MATMUL(invJ,dPsi))*r*wTria(l)/detJ
|
|
|
|
END DO
|
|
localK = localK*PI2
|
|
|
|
END FUNCTION elemKTria
|
|
|
|
!Compute 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):: detJ, f
|
|
REAL(8):: r
|
|
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
|
|
|
|
!Check if Xi is inside the element
|
|
PURE FUNCTION insideTria(Xi) RESULT(ins)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
LOGICAL:: ins
|
|
|
|
ins = Xi(1) >= 0.D0 .AND. &
|
|
Xi(2) >= 0.D0 .AND. &
|
|
1.D0 - Xi(1) - Xi(2) >= 0.D0
|
|
|
|
END FUNCTION insideTria
|
|
|
|
!Transform physical coordinates to element coordinates
|
|
PURE FUNCTION phy2logTria(self,r) RESULT(Xi)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylTria), INTENT(in):: self
|
|
REAL(8), INTENT(in):: r(1:3)
|
|
REAL(8):: Xi(1:3)
|
|
REAL(8):: deltaR(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)
|
|
pDer = self%partialDer(3, dPsi)
|
|
invJ = self%invJac(pDer)
|
|
detJ = self%detJac(pDer)
|
|
Xi = MATMUL(invJ,deltaR)/detJ
|
|
|
|
END FUNCTION phy2logTria
|
|
|
|
!Get the neighbour cell for a logical position Xi
|
|
SUBROUTINE neighbourElementTria(self, Xi, neighbourElement)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylTria), INTENT(in):: self
|
|
REAL(8), INTENT(in):: Xi(1:3)
|
|
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(neighbourElement)
|
|
SELECT CASE (nextInt)
|
|
CASE (1)
|
|
neighbourElement => self%e1
|
|
CASE (2)
|
|
neighbourElement => self%e2
|
|
CASE (3)
|
|
neighbourElement => self%e3
|
|
END SELECT
|
|
|
|
END SUBROUTINE neighbourElementTria
|
|
|
|
!Calculate volume for triangular element
|
|
PURE SUBROUTINE volumeTria(self)
|
|
USE moduleConstParam, ONLY: PI
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylTria), 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)
|
|
REAL(8):: r
|
|
|
|
self%volume = 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, 3)
|
|
r = DOT_PRODUCT(fPsi, self%r)
|
|
!Computes total volume of the cell
|
|
self%volume = r*detJ*PI !2PI*1/2
|
|
!Computes volume per node
|
|
self%n1%v = self%n1%v + fPsi(1)*self%volume
|
|
self%n2%v = self%n2%v + fPsi(2)*self%volume
|
|
self%n3%v = self%n3%v + fPsi(3)*self%volume
|
|
|
|
END SUBROUTINE volumeTria
|
|
|
|
!COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS
|
|
!Compute element Jacobian determinant
|
|
PURE FUNCTION detJ2DCyl(pDer) RESULT(dJ)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: pDer(1:3, 1:3)
|
|
REAL(8):: dJ
|
|
|
|
dJ = pDer(1,1)*pDer(2,2)-pDer(1,2)*pDer(2,1)
|
|
|
|
END FUNCTION detJ2DCyl
|
|
|
|
!Compute element Jacobian inverse matrix (without determinant)
|
|
PURE FUNCTION invJ2DCyl(pDer) RESULT(invJ)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: pDer(1:3, 1:3)
|
|
REAL(8):: invJ(1:3,1:3)
|
|
|
|
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
|
|
|
|
END FUNCTION invJ2DCyl
|
|
|
|
SUBROUTINE connectMesh2DCyl(self)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshGeneric), INTENT(inout):: self
|
|
INTEGER:: e, et
|
|
|
|
DO e = 1, self%numCells
|
|
!Connect Cell-Cell
|
|
DO et = 1, self%numCells
|
|
IF (e /= et) THEN
|
|
CALL connectCellCell(self%cells(e)%obj, self%cells(et)%obj)
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
SELECT TYPE(self)
|
|
TYPE IS(meshParticles)
|
|
!Connect Cell-Edge
|
|
DO et = 1, self%numEdges
|
|
CALL connectCellEdge(self%cells(e)%obj, self%edges(et)%obj)
|
|
|
|
END DO
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
END SUBROUTINE connectMesh2DCyl
|
|
|
|
!Selects type of elements to build connection
|
|
SUBROUTINE connectCellCell(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell), INTENT(inout):: elemA
|
|
CLASS(meshCell), INTENT(inout):: elemB
|
|
|
|
SELECT TYPE(elemA)
|
|
TYPE IS(meshCell2DCylQuad)
|
|
!Element A is a quadrilateral
|
|
SELECT TYPE(elemB)
|
|
TYPE IS(meshCell2DCylQuad)
|
|
!Element B is a quadrilateral
|
|
CALL connectQuadQuad(elemA, elemB)
|
|
|
|
TYPE IS(meshCell2DCylTria)
|
|
!Element B is a triangle
|
|
CALL connectQuadTria(elemA, elemB)
|
|
|
|
END SELECT
|
|
|
|
TYPE IS(meshCell2DCylTria)
|
|
!Element A is a Triangle
|
|
SELECT TYPE(elemB)
|
|
TYPE IS(meshCell2DCylQuad)
|
|
!Element B is a quadrilateral
|
|
CALL connectQuadTria(elemB, elemA)
|
|
|
|
TYPE IS(meshCell2DCylTria)
|
|
!Element B is a triangle
|
|
CALL connectTriaTria(elemA, elemB)
|
|
|
|
END SELECT
|
|
|
|
END SELECT
|
|
|
|
END SUBROUTINE connectCellCell
|
|
|
|
SUBROUTINE connectCellEdge(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell), INTENT(inout):: elemA
|
|
CLASS(meshEdge), INTENT(inout):: elemB
|
|
|
|
SELECT TYPE(elemB)
|
|
CLASS IS(meshEdge2DCyl)
|
|
SELECT TYPE(elemA)
|
|
TYPE IS(meshCell2DCylQuad)
|
|
!Element A is a quadrilateral
|
|
CALL connectQuadEdge(elemA, elemB)
|
|
|
|
TYPE IS(meshCell2DCylTria)
|
|
!Element A is a triangle
|
|
CALL connectTriaEdge(elemA, elemB)
|
|
|
|
END SELECT
|
|
|
|
END SELECT
|
|
|
|
END SUBROUTINE connectCellEdge
|
|
|
|
SUBROUTINE connectQuadQuad(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemA
|
|
CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemB
|
|
|
|
!Check direction 1
|
|
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
|
|
IF (elemA%n1%n == elemB%n4%n .AND. &
|
|
elemA%n2%n == elemB%n3%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
|
|
elemA%n2%n == elemB%n2%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
|
|
elemA%n2%n == elemB%n1%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n1%n == elemB%n1%n .AND. &
|
|
elemA%n2%n == elemB%n4%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e4 => elemA
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 2
|
|
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
|
|
IF (elemA%n2%n == elemB%n1%n .AND. &
|
|
elemA%n3%n == elemB%n4%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e4 => elemA
|
|
|
|
ELSEIF (elemA%n2%n == elemB%n4%n .AND. &
|
|
elemA%n3%n == elemB%n3%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
|
|
elemA%n3%n == elemB%n2%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
|
|
elemA%n3%n == elemB%n1%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 3
|
|
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
|
|
IF (elemA%n3%n == elemB%n2%n .AND. &
|
|
elemA%n4%n == elemB%n1%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n3%n == elemB%n1%n .AND. &
|
|
elemA%n4%n == elemB%n4%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e4 => elemA
|
|
|
|
ELSEIF (elemA%n3%n == elemB%n4%n .AND. &
|
|
elemA%n4%n == elemB%n3%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
|
|
elemA%n4%n == elemB%n2%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 4
|
|
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
|
|
IF (elemA%n4%n == elemB%n3%n .AND. &
|
|
elemA%n1%n == elemB%n2%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
|
|
elemA%n1%n == elemB%n1%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n4%n == elemB%n1%n .AND. &
|
|
elemA%n1%n == elemB%n4%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e4 => elemA
|
|
|
|
ELSEIF (elemA%n4%n == elemB%n4%n .AND. &
|
|
elemA%n1%n == elemB%n3%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE connectQuadQuad
|
|
|
|
SUBROUTINE connectQuadTria(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemA
|
|
CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemB
|
|
|
|
!Check direction 1
|
|
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
|
|
IF (elemA%n1%n == elemB%n1%n .AND. &
|
|
elemA%n2%n == elemB%n3%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
|
|
elemA%n2%n == elemB%n2%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
|
|
elemA%n2%n == elemB%n1%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 2
|
|
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
|
|
IF (elemA%n2%n == elemB%n1%n .AND. &
|
|
elemA%n3%n == elemB%n3%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
|
|
elemA%n3%n == elemB%n2%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
|
|
elemA%n3%n == elemB%n1%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 3
|
|
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
|
|
IF (elemA%n3%n == elemB%n1%n .AND. &
|
|
elemA%n4%n == elemB%n3%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
|
|
elemA%n4%n == elemB%n2%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
|
|
elemA%n4%n == elemB%n1%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 4
|
|
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
|
|
IF (elemA%n4%n == elemB%n1%n .AND. &
|
|
elemA%n1%n == elemB%n3%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n4%n == elemB%n3%n .AND. &
|
|
elemA%n1%n == elemB%n2%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
|
|
elemA%n1%n == elemB%n1%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE connectQuadTria
|
|
|
|
SUBROUTINE connectTriaTria(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemA
|
|
CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemB
|
|
|
|
!Check direction 1
|
|
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
|
|
IF (elemA%n1%n == elemB%n1%n .AND. &
|
|
elemA%n2%n == elemB%n3%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
|
|
elemA%n2%n == elemB%n1%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n1%n == elemB%n3%n .AND. &
|
|
elemA%n2%n == elemB%n2%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 2
|
|
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
|
|
IF (elemA%n2%n == elemB%n1%n .AND. &
|
|
elemA%n3%n == elemB%n3%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
|
|
elemA%n3%n == elemB%n1%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n2%n == elemB%n3%n .AND. &
|
|
elemA%n3%n == elemB%n2%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
END IF
|
|
|
|
|
|
END IF
|
|
|
|
!Check direction 3
|
|
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
|
|
IF (elemA%n3%n == elemB%n1%n .AND. &
|
|
elemA%n1%n == elemB%n3%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
|
|
elemA%n1%n == elemB%n1%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n3%n == elemB%n3%n .AND. &
|
|
elemA%n1%n == elemB%n2%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
END IF
|
|
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE connectTriaTria
|
|
|
|
SUBROUTINE connectQuadEdge(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylQuad), INTENT(inout), TARGET:: elemA
|
|
CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB
|
|
|
|
!Check direction 1
|
|
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
|
|
IF (elemA%n1%n == elemB%n1%n .AND. &
|
|
elemA%n2%n == elemB%n2%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
|
|
elemA%n2%n == elemB%n1%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
!Revers the normal to point inside the domain
|
|
elemB%normal = - elemB%normal
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 2
|
|
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
|
|
IF (elemA%n2%n == elemB%n1%n .AND. &
|
|
elemA%n3%n == elemB%n2%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
|
|
elemA%n3%n == elemB%n1%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
!Revers the normal to point inside the domain
|
|
elemB%normal = - elemB%normal
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 3
|
|
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
|
|
IF (elemA%n3%n == elemB%n1%n .AND. &
|
|
elemA%n4%n == elemB%n2%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
|
|
elemA%n4%n == elemB%n1%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
!Revers the normal to point inside the domain
|
|
elemB%normal = - elemB%normal
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 4
|
|
IF (.NOT. ASSOCIATED(elemA%e4)) THEN
|
|
IF (elemA%n4%n == elemB%n1%n .AND. &
|
|
elemA%n1%n == elemB%n2%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n4%n == elemB%n2%n .AND. &
|
|
elemA%n1%n == elemB%n1%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
!Revers the normal to point inside the domain
|
|
elemB%normal = - elemB%normal
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE connectQuadEdge
|
|
|
|
SUBROUTINE connectTriaEdge(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshCell2DCylTria), INTENT(inout), TARGET:: elemA
|
|
CLASS(meshEdge2DCyl), INTENT(inout), TARGET:: elemB
|
|
|
|
!Check direction 1
|
|
IF (.NOT. ASSOCIATED(elemA%e1)) THEN
|
|
IF (elemA%n1%n == elemB%n1%n .AND. &
|
|
elemA%n2%n == elemB%n2%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n1%n == elemB%n2%n .AND. &
|
|
elemA%n2%n == elemB%n1%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
!Revers the normal to point inside the domain
|
|
elemB%normal = - elemB%normal
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 2
|
|
IF (.NOT. ASSOCIATED(elemA%e2)) THEN
|
|
IF (elemA%n2%n == elemB%n1%n .AND. &
|
|
elemA%n3%n == elemB%n2%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n2%n == elemB%n2%n .AND. &
|
|
elemA%n3%n == elemB%n1%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
!Revers the normal to point inside the domain
|
|
elemB%normal = - elemB%normal
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
!Check direction 3
|
|
IF (.NOT. ASSOCIATED(elemA%e3)) THEN
|
|
IF (elemA%n3%n == elemB%n1%n .AND. &
|
|
elemA%n1%n == elemB%n2%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
|
|
elemA%n1%n == elemB%n1%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
!Revers the normal to point inside the domain
|
|
elemB%normal = - elemB%normal
|
|
|
|
END IF
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE connectTriaEdge
|
|
|
|
END MODULE moduleMesh2DCyl
|