JGonzalez/fpakc
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions
fpakc/src/modules/mesh/2DCyl/moduleMesh2DCylRead.f90
Jorge Gonzalez db6b0a2c03 Fixed an issue with reflection of particles in all geometries and also 
assigning the normal vector in 2D and 3D.

3D Cartesian geometry is working properly, although it needs testing.

Still issue with ionization boundary.
2021-03-22 12:39:34 +01:00

643 lines
16 KiB
Fortran
Raw Blame History

MODULE moduleMesh2DCylRead
USE moduleMesh
USE moduleMesh2DCyl
TYPE, EXTENDS(meshGeneric):: mesh2DCylGeneric
CONTAINS
PROCEDURE, PASS:: init => init2DCylMesh
PROCEDURE, PASS:: readMesh => readMesh2DCylGmsh
END TYPE
INTERFACE connected
MODULE PROCEDURE connectedVolVol, connectedVolEdge
END INTERFACE connected
CONTAINS
!Init mesh
SUBROUTINE init2DCylMesh(self, meshFormat)
USE moduleMesh
USE moduleErrors
IMPLICIT NONE
CLASS(mesh2DCylGeneric), INTENT(out):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: meshFormat
SELECT CASE(meshFormat)
CASE ("gmsh")
self%printOutput => printOutputGmsh
self%printColl => printCollGmsh
self%printEM => printEMGmsh
CASE DEFAULT
CALL criticalError("Mesh type " // meshFormat // " not supported.", "init2DCylMesh")
END SELECT
END SUBROUTINE init2DCylMesh
!Read mesh from gmsh file
SUBROUTINE readMesh2DCylGmsh(self, filename)
USE moduleBoundary
IMPLICIT NONE
CLASS(mesh2DCylGeneric), INTENT(inout):: self
CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
REAL(8):: r, z
INTEGER:: p(1:4)
INTEGER:: e=0, et=0, n=0, eTemp=0, elemType=0, bt = 0
INTEGER:: totalNumElem
INTEGER:: boundaryType
!Read msh
OPEN(10, FILE=TRIM(filename))
!Skip header
READ(10, *)
READ(10, *)
READ(10, *)
READ(10, *)
!Read number of nodes
READ(10, *) self%numNodes
!Allocate required matrices and vectors
ALLOCATE(self%nodes(1:self%numNodes))
ALLOCATE(self%K(1:self%numNodes,1:self%numNodes))
ALLOCATE(self%IPIV(1:self%numNodes,1:self%numNodes))
self%K = 0.D0
self%IPIV = 0
!Read nodes cartesian coordinates (x=z, y=r, z=null)
DO e=1, self%numNodes
READ(10, *) n, z, r
ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj)
CALL self%nodes(n)%obj%init(n, (/z, r, 0.D0 /))
END DO
!Skips comments
READ(10, *)
READ(10, *)
!Reads Totalnumber of elements
READ(10, *) TotalnumElem
!counts edges and volume elements
self%numEdges = 0
DO e=1, TotalnumElem
READ(10,*) eTemp, elemType
IF (elemType==1) THEN
self%numEdges=e
END IF
END DO
!Substract the number of edges to the total number of elements
!to obtain the number of volume elements
self%numVols = TotalnumElem - self%numEdges
!Allocates arrays
ALLOCATE(self%edges(1:self%numEdges))
ALLOCATE(self%vols(1:self%numVols))
!Go back to the beggining to read elements
DO e=1, totalNumElem
BACKSPACE(10)
END DO
!Reads edges
DO e=1, self%numEdges
READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
!Associate boundary condition procedure.
bt = getBoundaryId(boundaryType)
ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
CALL self%edges(e)%obj%init(n, p(1:2), bt, boundaryType)
END DO
!Read and initialize volumes
DO e=1, self%numVols
READ(10,*) n, elemType
BACKSPACE(10)
SELECT CASE(elemType)
CASE (2)
!Triangular element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:3)
ALLOCATE(meshVol2DCylTria:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:3))
CASE (3)
!Quadrilateral element
READ(10,*) n, elemType, eTemp, eTemp, eTemp, p(1:4)
ALLOCATE(meshVol2DCylQuad:: self%vols(e)%obj)
CALL self%vols(e)%obj%init(n - self%numEdges, p(1:4))
END SELECT
END DO
CLOSE(10)
!Build connectivity between elements
DO e = 1, self%numVols
!Connectivity between volumes
DO et = 1, self%numVols
IF (e /= et) THEN
CALL connected(self%vols(e)%obj, self%vols(et)%obj)
END IF
END DO
!Connectivity between vols and edges
DO et = 1, self%numEdges
CALL connected(self%vols(e)%obj, self%edges(et)%obj)
END DO
!Constructs the global K matrix
CALL constructGlobalK(self%K, self%vols(e)%obj)
END DO
END SUBROUTINE readMesh2DCylGmsh
!Selects type of elements to build connection
SUBROUTINE connectedVolVol(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 connectedQuadQuad(elemA, elemB)
TYPE IS(meshVol2DCylTria)
!Element B is a triangle
CALL connectedQuadTria(elemA, elemB)
END SELECT
TYPE IS(meshVol2DCylTria)
!Element A is a Triangle
SELECT TYPE(elemB)
TYPE IS(meshVol2DCylQuad)
!Element B is a quadrilateral
CALL connectedQuadTria(elemB, elemA)
TYPE IS(meshVol2DCylTria)
!Element B is a triangle
CALL connectedTriaTria(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolVol
SUBROUTINE connectedVolEdge(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 connectedQuadEdge(elemA, elemB)
TYPE IS(meshVol2DCylTria)
!Element A is a triangle
CALL connectedTriaEdge(elemA, elemB)
END SELECT
END SELECT
END SUBROUTINE connectedVolEdge
PURE FUNCTION coincidentNodes(nodesA, nodesB) RESULT(coincident)
IMPLICIT NONE
INTEGER, DIMENSION(1:2), INTENT(in):: nodesA, nodesB
LOGICAL:: coincident
INTEGER:: i
coincident = .FALSE.
DO i = 1, 2
IF (ANY(nodesA(i) == nodesB)) THEN
coincident = .TRUE.
ELSE
coincident = .FALSE.
EXIT
END IF
END DO
END FUNCTION coincidentNodes
SUBROUTINE connectedQuadQuad(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 connectedQuadQuad
SUBROUTINE connectedQuadTria(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 connectedQuadTria
SUBROUTINE connectedTriaTria(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 connectedTriaTria
SUBROUTINE connectedQuadEdge(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 connectedQuadEdge
SUBROUTINE connectedTriaEdge(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 connectedTriaEdge
SUBROUTINE constructGlobalK(K, elem)
IMPLICIT NONE
REAL(8), INTENT(inout):: K(1:,1:)
CLASS(meshVol), INTENT(in):: elem
REAL(8), ALLOCATABLE:: localK(:,:)
INTEGER:: nNodes, i, j
INTEGER, ALLOCATABLE:: n(:)
SELECT TYPE(elem)
TYPE IS(meshVol2DCylQuad)
nNodes = 4
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
ALLOCATE(n(1:nNodes))
n = (/ elem%n1%n, elem%n2%n, &
elem%n3%n, elem%n4%n /)
TYPE IS(meshVol2DCylTria)
nNodes = 3
ALLOCATE(localK(1:nNodes,1:nNodes))
localK = elem%elemK()
ALLOCATE(n(1:nNodes))
n = (/ elem%n1%n, elem%n2%n, elem%n3%n /)
CLASS DEFAULT
nNodes = 0
ALLOCATE(localK(1:1, 1:1))
localK = 0.D0
ALLOCATE(n(1:1))
n = 0
END SELECT
DO i = 1, nNodes
DO j = 1, nNodes
K(n(i), n(j)) = K(n(i), n(j)) + localK(i, j)
END DO
END DO
END SUBROUTINE constructGlobalK
END MODULE moduleMesh2DCylRead
Reference in a new issue View git blame Copy permalink