Jorge Gonzalez
e25b567d36
easy to use a file from a previous run without processing it into a plain text file. Although the previous method presented some updates for 1D small cases, this is quite easy to do with any type of simulations and, in the future, with different mesh formats.
1496 lines
41 KiB
Fortran
1496 lines
41 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
|
|
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
|
|
PROCEDURE, PASS:: init => initEdge2DCyl
|
|
PROCEDURE, PASS:: getNodes => getNodes2DCyl
|
|
PROCEDURE, PASS:: intersection => intersection2DCylEdge
|
|
PROCEDURE, PASS:: randPos => randPosEdge
|
|
|
|
END TYPE meshEdge2DCyl
|
|
|
|
TYPE, PUBLIC, ABSTRACT, EXTENDS(meshVol):: meshVol2DCyl
|
|
CONTAINS
|
|
PROCEDURE, PASS:: detJac => detJ2DCyl
|
|
PROCEDURE, PASS:: invJac => invJ2DCyl
|
|
PROCEDURE(dPsi_interface), DEFERRED, NOPASS:: dPsi
|
|
PROCEDURE(partialDer_interface), DEFERRED, PASS:: partialDer
|
|
|
|
END TYPE meshVol2DCyl
|
|
|
|
ABSTRACT INTERFACE
|
|
PURE FUNCTION dPsi_interface(xi) RESULT(dPsi)
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8), ALLOCATABLE:: dPsi(:,:)
|
|
|
|
END FUNCTION dPsi_interface
|
|
|
|
PURE SUBROUTINE partialDer_interface(self, dPsi, dz, dr)
|
|
IMPORT meshVol2DCyl
|
|
CLASS(meshVol2DCyl), INTENT(in):: self
|
|
REAL(8), INTENT(in):: dPsi(1:,1:)
|
|
REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
|
|
|
|
END SUBROUTINE partialDer_interface
|
|
|
|
END INTERFACE
|
|
|
|
!Quadrilateral volume element
|
|
TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylQuad
|
|
!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()
|
|
REAL(8):: arNodes(1:4) = 0.D0
|
|
|
|
CONTAINS
|
|
PROCEDURE, PASS:: init => initVolQuad2DCyl
|
|
PROCEDURE, PASS:: randPos => randPosVolQuad
|
|
PROCEDURE, PASS:: area => areaQuad
|
|
PROCEDURE, NOPASS:: fPsi => fPsiQuad
|
|
PROCEDURE, NOPASS:: dPsi => dPsiQuad
|
|
PROCEDURE, NOPASS:: dPsiXi1 => dPsiQuadXi1
|
|
PROCEDURE, NOPASS:: dPsiXi2 => dPsiQuadXi2
|
|
PROCEDURE, PASS:: partialDer => partialDerQuad
|
|
PROCEDURE, PASS:: elemK => elemKQuad
|
|
PROCEDURE, PASS:: elemF => elemFQuad
|
|
PROCEDURE, NOPASS:: weight => weightQuad
|
|
PROCEDURE, NOPASS:: inside => insideQuad
|
|
PROCEDURE, PASS:: scatter => scatterQuad
|
|
PROCEDURE, PASS:: gatherEF => gatherEFQuad
|
|
PROCEDURE, PASS:: getNodes => getNodesQuad
|
|
PROCEDURE, PASS:: phy2log => phy2logQuad
|
|
PROCEDURE, PASS:: nextElement => nextElementQuad
|
|
|
|
END TYPE meshVol2DCylQuad
|
|
|
|
!Triangular volume element
|
|
TYPE, PUBLIC, EXTENDS(meshVol2DCyl):: meshVol2DCylTria
|
|
!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()
|
|
REAL(8):: arNodes(1:3) = 0.D0
|
|
|
|
CONTAINS
|
|
PROCEDURE, PASS:: init => initVolTria2DCyl
|
|
PROCEDURE, PASS:: randPos => randPosVolTria
|
|
PROCEDURE, PASS:: area => areaTria
|
|
PROCEDURE, NOPASS:: fPsi => fPsiTria
|
|
PROCEDURE, NOPASS:: dPsi => dPsiTria
|
|
PROCEDURE, NOPASS:: dPsiXi1 => dPsiTriaXi1
|
|
PROCEDURE, NOPASS:: dPsiXi2 => dPsiTriaXi2
|
|
PROCEDURE, PASS:: partialDer => partialDerTria
|
|
PROCEDURE, PASS:: elemK => elemKTria
|
|
PROCEDURE, PASS:: elemF => elemFTria
|
|
PROCEDURE, NOPASS:: weight => weightTria
|
|
PROCEDURE, NOPASS:: inside => insideTria
|
|
PROCEDURE, PASS:: scatter => scatterTria
|
|
PROCEDURE, PASS:: gatherEF => gatherEFTria
|
|
PROCEDURE, PASS:: getNodes => getNodesTria
|
|
PROCEDURE, PASS:: phy2log => phy2logTria
|
|
PROCEDURE, PASS:: nextElement => nextElementTria
|
|
|
|
END TYPE meshVol2DCylTria
|
|
|
|
CONTAINS
|
|
!NODE FUNCTIONS
|
|
!Inits node element
|
|
SUBROUTINE initNode2DCyl(self, n, r)
|
|
USE moduleSpecies
|
|
USE moduleRefParam
|
|
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))
|
|
|
|
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
|
|
!Inits 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%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) RESULT(n)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshEdge2DCyl), INTENT(in):: self
|
|
INTEGER, ALLOCATABLE:: n(:)
|
|
|
|
ALLOCATE(n(1:2))
|
|
n = (/self%n1%n, self%n2%n /)
|
|
|
|
END FUNCTION getNodes2DCyl
|
|
|
|
PURE FUNCTION intersection2DCylEdge(self, r0) RESULT(r)
|
|
USE moduleMath
|
|
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
|
|
|
|
!Calculates a random position in edge
|
|
FUNCTION randPosEdge(self) RESULT(r)
|
|
USE moduleRandom
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshEdge2DCyl), INTENT(in):: self
|
|
REAL(8):: rnd
|
|
REAL(8):: dr, dz
|
|
REAL(8):: r(1:3)
|
|
|
|
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
|
|
!Inits quadrilateral element
|
|
SUBROUTINE initVolQuad2DCyl(self, n, p, nodes)
|
|
USE moduleRefParam
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), 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
|
|
|
|
self%n = n
|
|
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%area()
|
|
self%n1%v = self%n1%v + self%arNodes(1)
|
|
self%n2%v = self%n2%v + self%arNodes(2)
|
|
self%n3%v = self%n3%v + self%arNodes(3)
|
|
self%n4%v = self%n4%v + self%arNodes(4)
|
|
|
|
self%sigmaVrelMax = sigma_ref/L_ref**2
|
|
|
|
CALL OMP_INIT_LOCK(self%lock)
|
|
|
|
END SUBROUTINE initVolQuad2DCyl
|
|
|
|
!Computes element area
|
|
PURE SUBROUTINE areaQuad(self)
|
|
USE moduleConstParam, ONLY: PI8
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(inout):: self
|
|
REAL(8):: r, 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)*PI8 !4*2*pi
|
|
fPsi = self%fPsi(xi)
|
|
r = DOT_PRODUCT(fPsi,self%r)
|
|
self%volume = r*detJ
|
|
self%arNodes = fPsi*r*detJ
|
|
|
|
END SUBROUTINE areaQuad
|
|
|
|
!Computes element functions in point xi
|
|
PURE FUNCTION fPsiQuad(xi) RESULT(fPsi)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8),INTENT(in):: xi(1:3)
|
|
REAL(8), ALLOCATABLE:: fPsi(:)
|
|
|
|
ALLOCATE(fPsi(1:4))
|
|
|
|
fPsi(1) = (1.D0-xi(1))*(1.D0-xi(2))
|
|
fPsi(2) = (1.D0+xi(1))*(1.D0-xi(2))
|
|
fPsi(3) = (1.D0+xi(1))*(1.D0+xi(2))
|
|
fPsi(4) = (1.D0-xi(1))*(1.D0+xi(2))
|
|
fPsi = fPsi*0.25D0
|
|
|
|
END FUNCTION fPsiQuad
|
|
|
|
!Derivative element function at coordinates xi
|
|
PURE FUNCTION dPsiQuad(xi) RESULT(dPsi)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8), ALLOCATABLE:: dPsi(:,:)
|
|
|
|
ALLOCATE(dPsi(1:2,1:4))
|
|
|
|
dPsi(1,:) = dPsiQuadXi1(xi(2))
|
|
dPsi(2,:) = dPsiQuadXi2(xi(1))
|
|
|
|
END FUNCTION dPsiQuad
|
|
|
|
!Derivative element function (xi1)
|
|
PURE FUNCTION dPsiQuadXi1(xi2) RESULT(dPsiXi1)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8),INTENT(in):: xi2
|
|
REAL(8):: dPsiXi1(1:4)
|
|
|
|
dPsiXi1(1) = -(1.D0-xi2)
|
|
dPsiXi1(2) = (1.D0-xi2)
|
|
dPsiXi1(3) = (1.D0+xi2)
|
|
dPsiXi1(4) = -(1.D0+xi2)
|
|
dPsiXi1 = dPsiXi1*0.25D0
|
|
|
|
END FUNCTION dPsiQuadXi1
|
|
|
|
!Derivative element function (xi2)
|
|
PURE FUNCTION dPsiQuadXi2(xi1) RESULT(dPsiXi2)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8),INTENT(in):: xi1
|
|
REAL(8):: dPsiXi2(1:4)
|
|
|
|
dPsiXi2(1) = -(1.D0-xi1)
|
|
dPsiXi2(2) = -(1.D0+xi1)
|
|
dPsiXi2(3) = (1.D0+xi1)
|
|
dPsiXi2(4) = (1.D0-xi1)
|
|
dPsiXi2 = dPsiXi2*0.25D0
|
|
|
|
END FUNCTION dPsiQuadXi2
|
|
|
|
!Partial derivative in global coordinates
|
|
PURE SUBROUTINE partialDerQuad(self, dPsi, dz, dr)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(in):: self
|
|
REAL(8), INTENT(in):: dPsi(1:,1:)
|
|
REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
|
|
|
|
dz(1) = DOT_PRODUCT(dPsi(1,:),self%z)
|
|
dz(2) = DOT_PRODUCT(dPsi(2,:),self%z)
|
|
dr(1) = DOT_PRODUCT(dPsi(1,:),self%r)
|
|
dr(2) = DOT_PRODUCT(dPsi(2,:),self%r)
|
|
|
|
END SUBROUTINE partialDerQuad
|
|
|
|
!Random position in quadrilateral volume
|
|
FUNCTION randPosVolQuad(self) RESULT(r)
|
|
USE moduleRandom
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
REAL(8):: xii(1:3)
|
|
REAL(8), ALLOCATABLE:: fPsi(:)
|
|
|
|
xii(1) = random(-1.D0, 1.D0)
|
|
xii(2) = random(-1.D0, 1.D0)
|
|
xii(3) = 0.D0
|
|
|
|
fPsi = self%fPsi(xii)
|
|
|
|
r(1) = DOT_PRODUCT(fPsi, self%z)
|
|
r(2) = DOT_PRODUCT(fPsi, self%r)
|
|
r(3) = 0.D0
|
|
|
|
END FUNCTION randposVolQuad
|
|
|
|
!Computes element local stiffness matrix
|
|
PURE FUNCTION elemKQuad(self) RESULT(localK)
|
|
USE moduleConstParam, ONLY: PI2
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(in):: self
|
|
REAL(8), ALLOCATABLE:: localK(:,:)
|
|
REAL(8):: r, xi(1:3)
|
|
REAL(8):: fPsi(1:4), dPsi(1:2,1:4)
|
|
REAL(8):: invJ(1:2,1:2), detJ
|
|
INTEGER:: l, m
|
|
|
|
ALLOCATE(localK(1:4, 1:4))
|
|
localK=0.D0
|
|
xi=0.D0
|
|
!Start 2D Gauss Quad Integral
|
|
DO l=1, 3
|
|
xi(2) = corQuad(l)
|
|
dPsi(1,:) = self%dPsiXi1(xi(2))
|
|
DO m = 1, 3
|
|
xi(1) = corQuad(m)
|
|
dPsi(2,:) = self%dPsiXi2(xi(1))
|
|
fPsi = self%fPsi(xi)
|
|
detJ = self%detJac(xi,dPsi)
|
|
invJ = self%invJac(xi,dPsi)
|
|
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, source) RESULT(localF)
|
|
USE moduleConstParam, ONLY: PI2
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(in):: self
|
|
REAL(8), INTENT(in):: source(1:)
|
|
REAL(8), ALLOCATABLE:: localF(:)
|
|
REAL(8):: r, xi(1:3)
|
|
REAL(8):: fPsi(1:4)
|
|
REAL(8):: detJ, f
|
|
INTEGER:: l, m
|
|
|
|
ALLOCATE(localF(1:4))
|
|
localF = 0.D0
|
|
xi = 0.D0
|
|
DO l=1, 3
|
|
xi(1) = corQuad(l)
|
|
DO m = 1, 3
|
|
xi(2) = corQuad(m)
|
|
detJ = self%detJac(xi)
|
|
fPsi = self%fPsi(xi)
|
|
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
|
|
|
|
!Computes weights in the element nodes
|
|
PURE FUNCTION weightQuad(xi) RESULT(w)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8):: w(1:4)
|
|
|
|
w = fPsiQuad(xi)
|
|
|
|
END FUNCTION weightQuad
|
|
|
|
!Checks if a particle is inside a quad 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
|
|
|
|
!Scatter properties of particle into element nodes
|
|
SUBROUTINE scatterQuad(self, part)
|
|
USE moduleMath
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(in):: self
|
|
CLASS(particle), INTENT(in):: part
|
|
TYPE(outputNode), POINTER:: vertex
|
|
REAL(8):: w_p(1:4)
|
|
REAL(8):: tensorS(1:3,1:3)
|
|
|
|
w_p = self%weight(part%xi)
|
|
tensorS = outerProduct(part%v, part%v)
|
|
|
|
vertex => self%n1%output(part%species%n)
|
|
vertex%den = vertex%den + part%weight*w_p(1)
|
|
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:)
|
|
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
|
|
|
|
vertex => self%n2%output(part%species%n)
|
|
vertex%den = vertex%den + part%weight*w_p(2)
|
|
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:)
|
|
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
|
|
|
|
vertex => self%n3%output(part%species%n)
|
|
vertex%den = vertex%den + part%weight*w_p(3)
|
|
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:)
|
|
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
|
|
|
|
vertex => self%n4%output(part%species%n)
|
|
vertex%den = vertex%den + part%weight*w_p(4)
|
|
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(4)*part%v(:)
|
|
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(4)*tensorS
|
|
|
|
END SUBROUTINE scatterQuad
|
|
|
|
!Gathers the electric field at position xi
|
|
PURE FUNCTION gatherEFQuad(self,xi) RESULT(EF)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(in):: self
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8):: dPsi(1:2,1:4)
|
|
REAL(8):: dPsiR(1:2,1:4)!Derivative of shpae functions in global coordinates
|
|
REAL(8):: invJ(1:2,1:2), detJ
|
|
REAL(8):: phi(1:4)
|
|
REAL(8):: EF(1:3)
|
|
|
|
phi = (/self%n1%emData%phi, &
|
|
self%n2%emData%phi, &
|
|
self%n3%emData%phi, &
|
|
self%n4%emData%phi /)
|
|
|
|
dPsi = self%dPsi(xi)
|
|
detJ = self%detJac(xi,dPsi)
|
|
invJ = self%invJac(xi,dPsi)
|
|
dPsiR = MATMUL(invJ, dPsi)/detJ
|
|
EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi)
|
|
EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi)
|
|
EF(3) = 0.D0
|
|
|
|
END FUNCTION gatherEFQuad
|
|
|
|
!Gets nodes from quadrilateral element
|
|
PURE FUNCTION getNodesQuad(self) RESULT(n)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(in):: self
|
|
INTEGER, ALLOCATABLE:: n(:)
|
|
|
|
ALLOCATE(n(1:4))
|
|
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(xN)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(in):: self
|
|
REAL(8), INTENT(in):: r(1:3)
|
|
REAL(8):: xN(1:3)
|
|
REAL(8):: xO(1:3), detJ, invJ(1:2,1:2), f(1:2)
|
|
REAL(8):: dPsi(1:2,1:4), fPsi(1:4)
|
|
REAL(8):: conv
|
|
|
|
!Iterative newton method to transform coordinates
|
|
conv=1.D0
|
|
xO=0.D0
|
|
|
|
DO WHILE(conv>1.D-4)
|
|
dPsi = self%dPsi(xO)
|
|
invJ = self%invJac(xO, dPsi)
|
|
fPsi = self%fPsi(xO)
|
|
f(1) = DOT_PRODUCT(fPsi,self%z)-r(1)
|
|
f(2) = DOT_PRODUCT(fPsi,self%r)-r(2)
|
|
detJ = self%detJac(xO,dPsi)
|
|
xN(1:2)=xO(1:2) - MATMUL(invJ, f)/detJ
|
|
conv=MAXVAL(DABS(xN-xO),1)
|
|
xO=xN
|
|
|
|
END DO
|
|
|
|
END FUNCTION phy2logQuad
|
|
|
|
!Gets the next element for a logical position xi
|
|
SUBROUTINE nextElementQuad(self, xi, nextElement)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(in):: self
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
CLASS(meshElement), POINTER, INTENT(out):: nextElement
|
|
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)
|
|
SELECT CASE (nextInt)
|
|
CASE (1)
|
|
nextElement => self%e1
|
|
CASE (2)
|
|
nextElement => self%e2
|
|
CASE (3)
|
|
nextElement => self%e3
|
|
CASE (4)
|
|
nextElement => self%e4
|
|
END SELECT
|
|
|
|
END SUBROUTINE nextElementQuad
|
|
|
|
!TRIA ELEMENT
|
|
!Init tria element
|
|
SUBROUTINE initVolTria2DCyl(self, n, p, nodes)
|
|
USE moduleRefParam
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), 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 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%area()
|
|
self%n1%v = self%n1%v + self%arNodes(1)
|
|
self%n2%v = self%n2%v + self%arNodes(2)
|
|
self%n3%v = self%n3%v + self%arNodes(3)
|
|
|
|
self%sigmaVrelMax = sigma_ref/L_ref**2
|
|
|
|
CALL OMP_INIT_LOCK(self%lock)
|
|
|
|
END SUBROUTINE initVolTria2DCyl
|
|
|
|
!Random position in quadrilateral volume
|
|
FUNCTION randPosVolTria(self) RESULT(r)
|
|
USE moduleRandom
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(in):: self
|
|
REAL(8):: r(1:3)
|
|
REAL(8):: xii(1:3)
|
|
REAL(8), ALLOCATABLE:: fPsi(:)
|
|
|
|
xii(1) = random( 0.D0, 1.D0)
|
|
xii(2) = random( 0.D0, 1.D0)
|
|
xii(3) = 0.D0
|
|
|
|
ALLOCATE(fPsi(1:3))
|
|
fPsi = self%fPsi(xii)
|
|
|
|
r(1) = DOT_PRODUCT(fPsi, self%z)
|
|
r(2) = DOT_PRODUCT(fPsi, self%r)
|
|
r(3) = 0.D0
|
|
|
|
END FUNCTION randposVolTria
|
|
|
|
!Calculates area for triangular element
|
|
PURE SUBROUTINE areaTria(self)
|
|
USE moduleConstParam, ONLY: PI
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(inout):: self
|
|
REAL(8):: r, 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)*PI !2PI*1/2
|
|
fPsi = self%fPsi(xi)
|
|
r = DOT_PRODUCT(fPsi,self%r)
|
|
self%volume = r*detJ
|
|
self%arNodes = fPsi*r*detJ
|
|
|
|
END SUBROUTINE areaTria
|
|
|
|
!Shape functions for triangular element
|
|
PURE FUNCTION fPsiTria(xi) RESULT(fPsi)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8), ALLOCATABLE:: fPsi(:)
|
|
|
|
ALLOCATE(fPsi(1:3))
|
|
|
|
fPsi(1) = 1.D0 - xi(1) - xi(2)
|
|
fPsi(2) = xi(1)
|
|
fPsi(3) = xi(2)
|
|
|
|
END FUNCTION fPsiTria
|
|
|
|
!Derivative element function at coordinates xi
|
|
PURE FUNCTION dPsiTria(xi) RESULT(dPsi)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8), ALLOCATABLE:: dPsi(:,:)
|
|
|
|
ALLOCATE(dPsi(1:2,1:3))
|
|
|
|
dPsi(1,:) = dPsiTriaXi1(xi(2))
|
|
dPsi(2,:) = dPsiTriaXi2(xi(1))
|
|
|
|
END FUNCTION dPsiTria
|
|
|
|
!Derivative element function (xi1)
|
|
PURE FUNCTION dPsiTriaXi1(xi2) RESULT(dPsiXi1)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: xi2
|
|
REAL(8):: dPsiXi1(1:3)
|
|
|
|
dPsiXi1(1) = -1.D0
|
|
dPsiXi1(2) = 1.D0
|
|
dPsiXi1(3) = 0.D0
|
|
|
|
END FUNCTION dPsiTriaXi1
|
|
|
|
!Derivative element function (xi1)
|
|
PURE FUNCTION dPsiTriaXi2(xi1) RESULT(dPsiXi2)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: xi1
|
|
REAL(8):: dPsiXi2(1:3)
|
|
|
|
dPsiXi2(1) = -1.D0
|
|
dPsiXi2(2) = 0.D0
|
|
dPsiXi2(3) = 1.D0
|
|
|
|
END FUNCTION dPsiTriaXi2
|
|
|
|
PURE SUBROUTINE partialDerTria(self, dPsi, dz, dr)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(in):: self
|
|
REAL(8), INTENT(in):: dPsi(1:,1:)
|
|
REAL(8), INTENT(out), DIMENSION(1:2):: dz, dr
|
|
|
|
dz(1) = DOT_PRODUCT(dPsi(1,:),self%z)
|
|
dz(2) = DOT_PRODUCT(dPsi(2,:),self%z)
|
|
dr(1) = DOT_PRODUCT(dPsi(1,:),self%r)
|
|
dr(2) = DOT_PRODUCT(dPsi(2,:),self%r)
|
|
|
|
END SUBROUTINE partialDerTria
|
|
|
|
!Computes element local stiffness matrix
|
|
PURE FUNCTION elemKTria(self) RESULT(localK)
|
|
USE moduleConstParam, ONLY: PI2
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(in):: self
|
|
REAL(8), ALLOCATABLE:: localK(:,:)
|
|
REAL(8):: r, xi(1:3)
|
|
REAL(8):: fPsi(1:3), dPsi(1:2,1:3)
|
|
REAL(8):: invJ(1:2,1:2), detJ
|
|
INTEGER:: l
|
|
|
|
ALLOCATE(localK(1:4, 1:4))
|
|
localK=0.D0
|
|
xi=0.D0
|
|
!Start 2D Gauss Quad Integral
|
|
DO l=1, 4
|
|
xi(1) = xi1Tria(l)
|
|
xi(2) = xi2Tria(l)
|
|
dPsi = self%dPsi(xi)
|
|
detJ = self%detJac(xi,dPsi)
|
|
invJ = self%invJac(xi,dPsi)
|
|
fPsi = self%fPsi(xi)
|
|
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, source) RESULT(localF)
|
|
USE moduleConstParam, ONLY: PI2
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(in):: self
|
|
REAL(8), INTENT(in):: source(1:)
|
|
REAL(8), ALLOCATABLE:: localF(:)
|
|
REAL(8):: fPsi(1:3)
|
|
REAL(8):: r, xi(1:3)
|
|
REAL(8):: detJ, f
|
|
INTEGER:: l
|
|
|
|
ALLOCATE(localF(1:3))
|
|
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)
|
|
fPsi = self%fPsi(xi)
|
|
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
|
|
|
|
!Computes weights in the element nodes
|
|
PURE FUNCTION weightTria(xi) RESULT(w)
|
|
IMPLICIT NONE
|
|
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8), ALLOCATABLE:: w(:)
|
|
|
|
w = fPsiTria(xi)
|
|
|
|
END FUNCTION weightTria
|
|
|
|
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
|
|
|
|
!Scatter properties of particles into element
|
|
SUBROUTINE scatterTria(self, part)
|
|
USE moduleMath
|
|
USE moduleSpecies
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(in):: self
|
|
CLASS(particle), INTENT(in):: part
|
|
TYPE(outputNode), POINTER:: vertex
|
|
REAL(8):: w_p(1:3)
|
|
REAL(8):: tensorS(1:3,1:3)
|
|
|
|
w_p = self%weight(part%xi)
|
|
tensorS = outerProduct(part%v, part%v)
|
|
|
|
vertex => self%n1%output(part%species%n)
|
|
vertex%den = vertex%den + part%weight*w_p(1)
|
|
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(1)*part%v(:)
|
|
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(1)*tensorS
|
|
|
|
vertex => self%n2%output(part%species%n)
|
|
vertex%den = vertex%den + part%weight*w_p(2)
|
|
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(2)*part%v(:)
|
|
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(2)*tensorS
|
|
|
|
vertex => self%n3%output(part%species%n)
|
|
vertex%den = vertex%den + part%weight*w_p(3)
|
|
vertex%mom(:) = vertex%mom(:) + part%weight*w_p(3)*part%v(:)
|
|
vertex%tensorS(:,:) = vertex%tensorS(:,:) + part%weight*w_p(3)*tensorS
|
|
|
|
END SUBROUTINE scatterTria
|
|
|
|
!Gathers the electric field at position xi
|
|
PURE FUNCTION gatherEFTria(self,xi) RESULT(EF)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(in):: self
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8):: dPsi(1:2,1:3)
|
|
REAL(8):: dPsiR(1:2,1:3)!Derivative of shpae functions in global coordinates
|
|
REAL(8):: invJ(1:2,1:2), detJ
|
|
REAL(8):: phi(1:3)
|
|
REAL(8):: EF(1:3)
|
|
|
|
phi = (/self%n1%emData%phi, &
|
|
self%n2%emData%phi, &
|
|
self%n3%emData%phi /)
|
|
|
|
dPsi = self%dPsi(xi)
|
|
detJ = self%detJac(xi,dPsi)
|
|
invJ = self%invJac(xi,dPsi)
|
|
dPsiR = MATMUL(invJ, dPsi)/detJ
|
|
EF(1) = -DOT_PRODUCT(dPsiR(1,:), phi)
|
|
EF(2) = -DOT_PRODUCT(dPsiR(2,:), phi)
|
|
EF(3) = 0.D0
|
|
|
|
END FUNCTION gatherEFTria
|
|
|
|
!Gets node indexes from triangular element
|
|
PURE FUNCTION getNodesTria(self) RESULT(n)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(in):: self
|
|
INTEGER, ALLOCATABLE:: n(:)
|
|
|
|
ALLOCATE(n(1:3))
|
|
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
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(in):: self
|
|
REAL(8), INTENT(in):: r(1:3)
|
|
REAL(8):: xi(1:3)
|
|
REAL(8):: invJ(1:2,1:2), detJ
|
|
REAL(8):: deltaR(1:2)
|
|
REAL(8):: dPsi(1:2,1:3)
|
|
|
|
!Direct method to convert coordinates
|
|
xi = 0.D0 !Irrelevant, required for input
|
|
deltaR = (/ r(1) - self%z(1), r(2) - self%r(1) /)
|
|
dPsi = self%dPsi(xi)
|
|
invJ = self%invJac(xi, dPsi)
|
|
detJ = self%detJac(xi, dPsi)
|
|
xi(1:2) = MATMUL(invJ,deltaR)/detJ
|
|
|
|
END FUNCTION phy2logTria
|
|
|
|
SUBROUTINE nextElementTria(self, xi, nextElement)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylTria), INTENT(in):: self
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
CLASS(meshElement), POINTER, INTENT(out):: nextElement
|
|
REAL(8):: xiArray(1:3)
|
|
INTEGER:: nextInt
|
|
|
|
xiArray = (/ xi(2), 1.D0-xi(1)-xi(2), xi(1) /)
|
|
nextInt = MINLOC(xiArray,1)
|
|
NULLIFY(nextElement)
|
|
SELECT CASE (nextInt)
|
|
CASE (1)
|
|
nextElement => self%e1
|
|
CASE (2)
|
|
nextElement => self%e2
|
|
CASE (3)
|
|
nextElement => self%e3
|
|
END SELECT
|
|
|
|
END SUBROUTINE nextElementTria
|
|
|
|
!COMMON FUNCTIONS FOR CYLINDRICAL VOLUME ELEMENTS
|
|
!Computes element Jacobian determinant
|
|
PURE FUNCTION detJ2DCyl(self, xi, dPsi_in) RESULT(dJ)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCyl), INTENT(in):: self
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
|
|
REAL(8), ALLOCATABLE:: dPsi(:,:)
|
|
REAL(8):: dJ
|
|
REAL(8):: dz(1:2), dr(1:2)
|
|
|
|
IF(PRESENT(dPsi_in)) THEN
|
|
dPsi = dPsi_in
|
|
|
|
ELSE
|
|
dPsi = self%dPsi(xi)
|
|
|
|
END IF
|
|
|
|
CALL self%partialDer(dPsi, dz, dr)
|
|
dJ = dz(1)*dr(2)-dz(2)*dr(1)
|
|
|
|
END FUNCTION detJ2DCyl
|
|
|
|
!Computes element Jacobian inverse matrix (without determinant)
|
|
PURE FUNCTION invJ2DCyl(self,xi,dPsi_in) RESULT(invJ)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCyl), INTENT(in):: self
|
|
REAL(8), INTENT(in):: xi(1:3)
|
|
REAL(8), INTENT(in), OPTIONAL:: dPsi_in(1:,1:)
|
|
REAL(8), ALLOCATABLE:: dPsi(:,:)
|
|
REAL(8):: dz(1:2), dr(1:2)
|
|
REAL(8):: invJ(1:2,1:2)
|
|
|
|
IF(PRESENT(dPsi_in)) THEN
|
|
dPsi=dPsi_in
|
|
|
|
ELSE
|
|
dPsi = self%dPsi(xi)
|
|
|
|
END IF
|
|
|
|
CALL self%partialDer(dPsi, dz, dr)
|
|
invJ(1,:) = (/ dr(2), -dz(2) /)
|
|
invJ(2,:) = (/ -dr(1), dz(1) /)
|
|
|
|
END FUNCTION invJ2DCyl
|
|
|
|
SUBROUTINE connectMesh2DCyl(self)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshGeneric), INTENT(inout):: self
|
|
INTEGER:: e, et
|
|
|
|
DO e = 1, self%numVols
|
|
!Connect Vol-Vol
|
|
DO et = 1, self%numVols
|
|
IF (e /= et) THEN
|
|
CALL connectVolVol(self%vols(e)%obj, self%vols(et)%obj)
|
|
|
|
END IF
|
|
|
|
END DO
|
|
|
|
SELECT TYPE(self)
|
|
TYPE IS(meshParticles)
|
|
!Connect Vol-Edge
|
|
DO et = 1, self%numEdges
|
|
CALL connectVolEdge(self%vols(e)%obj, self%edges(et)%obj)
|
|
|
|
END DO
|
|
|
|
END SELECT
|
|
|
|
END DO
|
|
|
|
END SUBROUTINE connectMesh2DCyl
|
|
|
|
!Selects type of elements to build connection
|
|
SUBROUTINE connectVolVol(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol), INTENT(inout):: elemA
|
|
CLASS(meshVol), INTENT(inout):: elemB
|
|
|
|
SELECT TYPE(elemA)
|
|
TYPE IS(meshVol2DCylQuad)
|
|
!Element A is a quadrilateral
|
|
SELECT TYPE(elemB)
|
|
TYPE IS(meshVol2DCylQuad)
|
|
!Element B is a quadrilateral
|
|
CALL connectQuadQuad(elemA, elemB)
|
|
|
|
TYPE IS(meshVol2DCylTria)
|
|
!Element B is a triangle
|
|
CALL connectQuadTria(elemA, elemB)
|
|
|
|
END SELECT
|
|
|
|
TYPE IS(meshVol2DCylTria)
|
|
!Element A is a Triangle
|
|
SELECT TYPE(elemB)
|
|
TYPE IS(meshVol2DCylQuad)
|
|
!Element B is a quadrilateral
|
|
CALL connectQuadTria(elemB, elemA)
|
|
|
|
TYPE IS(meshVol2DCylTria)
|
|
!Element B is a triangle
|
|
CALL connectTriaTria(elemA, elemB)
|
|
|
|
END SELECT
|
|
|
|
END SELECT
|
|
|
|
END SUBROUTINE connectVolVol
|
|
|
|
SUBROUTINE connectVolEdge(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol), INTENT(inout):: elemA
|
|
CLASS(meshEdge), INTENT(inout):: elemB
|
|
|
|
SELECT TYPE(elemB)
|
|
CLASS IS(meshEdge2DCyl)
|
|
SELECT TYPE(elemA)
|
|
TYPE IS(meshVol2DCylQuad)
|
|
!Element A is a quadrilateral
|
|
CALL connectQuadEdge(elemA, elemB)
|
|
|
|
TYPE IS(meshVol2DCylTria)
|
|
!Element A is a triangle
|
|
CALL connectTriaEdge(elemA, elemB)
|
|
|
|
END SELECT
|
|
|
|
END SELECT
|
|
|
|
END SUBROUTINE connectVolEdge
|
|
|
|
SUBROUTINE connectQuadQuad(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA
|
|
CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemB
|
|
|
|
!Check direction 1
|
|
IF (.NOT. ASSOCIATED(elemA%e1) .AND. &
|
|
elemA%n1%n == elemB%n4%n .AND. &
|
|
elemA%n2%n == elemB%n3%n) THEN
|
|
elemA%e1 => elemB
|
|
elemB%e3 => elemA
|
|
|
|
END IF
|
|
|
|
!Check direction 2
|
|
IF (.NOT. ASSOCIATED(elemA%e2) .AND. &
|
|
elemA%n2%n == elemB%n1%n .AND. &
|
|
elemA%n3%n == elemB%n4%n) THEN
|
|
elemA%e2 => elemB
|
|
elemB%e4 => elemA
|
|
|
|
END IF
|
|
|
|
!Check direction 3
|
|
IF (.NOT. ASSOCIATED(elemA%e3) .AND. &
|
|
elemA%n3%n == elemB%n2%n .AND. &
|
|
elemA%n4%n == elemB%n1%n) THEN
|
|
elemA%e3 => elemB
|
|
elemB%e1 => elemA
|
|
|
|
END IF
|
|
|
|
!Check direction 4
|
|
IF (.NOT. ASSOCIATED(elemA%e4) .AND. &
|
|
elemA%n4%n == elemB%n3%n .AND. &
|
|
elemA%n1%n == elemB%n2%n) THEN
|
|
elemA%e4 => elemB
|
|
elemB%e2 => elemA
|
|
|
|
END IF
|
|
|
|
END SUBROUTINE connectQuadQuad
|
|
|
|
SUBROUTINE connectQuadTria(elemA, elemB)
|
|
IMPLICIT NONE
|
|
|
|
CLASS(meshVol2DCylQuad), INTENT(inout), TARGET:: elemA
|
|
CLASS(meshVol2DCylTria), 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(meshVol2DCylTria), INTENT(inout), TARGET:: elemA
|
|
CLASS(meshVol2DCylTria), 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(meshVol2DCylQuad), 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(meshVol2DCylTria), 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
|