MODULE moduleMeshCylRead
  USE moduleMesh
  USE moduleMeshCyl
  USE moduleMeshCylBoundary

  TYPE, EXTENDS(meshGeneric):: meshCylGeneric
    CONTAINS
      PROCEDURE, PASS:: readMesh => readMeshCyl
      PROCEDURE, PASS:: printOutput => printOutputCyl
      PROCEDURE, PASS:: printColl   => printCollisionsCyl
      PROCEDURE, PASS:: printEM     => printEMCyl

  END TYPE

  INTERFACE connected
    MODULE PROCEDURE connectedVolVol, connectedVolEdge

  END INTERFACE connected

  CONTAINS
    SUBROUTINE readMeshCyl(self, filename)
      USE moduleRefParam
      USE moduleBoundary
      IMPLICIT NONE

      CLASS(meshCylGeneric), INTENT(out):: 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(meshNodeCyl:: self%nodes(n)%obj)
        CALL self%nodes(n)%obj%init(n, (/r, z, 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.
        !TODO: move to subroutine
        bt = getBoundaryId(boundaryType)
        SELECT CASE(boundary(bt)%obj%boundaryType)
        CASE ('reflection')
          ALLOCATE(meshEdgeCylRef:: self%edges(e)%obj)

        CASE ('absorption')
          ALLOCATE(meshEdgeCylAbs:: self%edges(e)%obj)

        CASE ('axis')
          ALLOCATE(meshEdgeCylAxis:: self%edges(e)%obj)

        END SELECT

        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(meshVolCylTria:: 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(meshVolCylQuad:: 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
        IF (e /= et) THEN
          DO et = 1, self%numVols
            CALL connected(self%vols(e)%obj, self%vols(et)%obj)

          END DO
        END IF

        !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

    !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(meshVolCylQuad)
        !Element A is a quadrilateral
        SELECT TYPE(elemB)
        TYPE IS(meshVolCylQuad)
          !Element B is a quadrilateral
          CALL connectedQuadQuad(elemA, elemB)

        TYPE IS(meshVolCylTria)
          !Element B is a triangle
          CALL connectedQuadTria(elemA, elemB)

        END SELECT

      TYPE IS(meshVolCylTria)
        !Element A is a Triangle
        SELECT TYPE(elemB)
        TYPE IS(meshVolCylQuad)
          !Element B is a quadrilateral
          CALL connectedQuadTria(elemB, elemA)

        TYPE IS(meshVolCylTria)
          !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(meshEdgeCyl)
        SELECT TYPE(elemA)
        TYPE IS(meshVolCylQuad)
          !Element A is a quadrilateral
          CALL connectedQuadEdge(elemA, elemB)

        TYPE IS(meshVolCylTria)
          !Element A is a triangle
          CALL connectedTriaEdge(elemA, elemB)

        END SELECT

      END SELECT

    END SUBROUTINE connectedVolEdge

    SUBROUTINE connectedQuadQuad(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshVolCylQuad), INTENT(inout), TARGET:: elemA
      CLASS(meshVolCylQuad), 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(meshVolCylQuad), INTENT(inout), TARGET:: elemA
      CLASS(meshVolCylTria), 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(meshVolCylTria), INTENT(inout), TARGET:: elemA
      CLASS(meshVolCylTria), 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(meshVolCylQuad), INTENT(inout), TARGET:: elemA
      CLASS(meshEdgeCyl),    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%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%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%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%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 connectedQuadEdge

    SUBROUTINE connectedTriaEdge(elemA, elemB)
      IMPLICIT NONE

      CLASS(meshVolCylTria), INTENT(inout), TARGET:: elemA
      CLASS(meshEdgeCyl),    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%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%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%n1%n == elemB%n2%n) THEN
          elemA%e3 => elemB
          elemB%e2 => elemA

        ELSEIF (elemA%n3%n == elemB%n2%n .AND. &
                elemA%n1%n == elemB%n1%n) THEN
          elemA%e3 => elemB
          elemB%e1 => elemA

        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(meshVolCylQuad)
        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(meshVolCylTria)
        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

    SUBROUTINE printOutputCyl(self, t)
      USE moduleRefParam
      USE moduleSpecies
      USE moduleOutput
      IMPLICIT NONE

      CLASS(meshCylGeneric), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: n, i
      TYPE(outputFormat):: output(1:self%numNodes)
      REAL(8):: time
      CHARACTER(:), ALLOCATABLE:: fileName
      CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations

      time = DBLE(t)*tau*ti_ref

      DO i = 1, nSpecies
        WRITE(tstring, '(I6.6)') t
        fileName='OUTPUT_' // tstring// '_' // species(i)%obj%name // '.msh'
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (60, file = path // folder // '/' //  fileName)
        WRITE(60, "(A)") '$MeshFormat'
        WRITE(60, "(A)") '2.2 0 8'
        WRITE(60, "(A)") '$EndMeshFormat'
        WRITE(60, "(A)") '$NodeData'
        WRITE(60, "(A)") '1'
        WRITE(60, "(A)") '"Density (m^-3)"'
        WRITE(60, *) 1
        WRITE(60, *) time
        WRITE(60, *) 3
        WRITE(60, *) t
        WRITE(60, *) 1
        WRITE(60, *) self%numNodes
        DO n=1, self%numNodes
          CALL calculateOutput(self%nodes(n)%obj%output(i), output(n), self%nodes(n)%obj%v, species(i)%obj)
          WRITE(60, "(I6,ES20.6E3)") n, output(n)%density
        END DO
        WRITE(60, "(A)") '$EndNodeData'
        WRITE(60, "(A)") '$NodeData'
        WRITE(60, "(A)") '1'
        WRITE(60, "(A)") '"Velocity (m/s)"'
        WRITE(60, *) 1
        WRITE(60, *) time
        WRITE(60, *) 3
        WRITE(60, *) t
        WRITE(60, *) 3
        WRITE(60, *) self%numNodes
        DO n=1, self%numNodes
          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%velocity
        END DO
        WRITE(60, "(A)") '$EndNodeData'
        WRITE(60, "(A)") '$NodeData'
        WRITE(60, "(A)") '1'
        WRITE(60, "(A)") '"Pressure (Pa)"'
        WRITE(60, *) 1
        WRITE(60, *) time
        WRITE(60, *) 3
        WRITE(60, *) t
        WRITE(60, *) 1
        WRITE(60, *) self%numNodes
        DO n=1, self%numNodes
          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%pressure
        END DO
        WRITE(60, "(A)") '$EndNodeData'
        WRITE(60, "(A)") '$NodeData'
        WRITE(60, "(A)") '1'
        WRITE(60, "(A)") '"Temperature (K)"'
        WRITE(60, *) 1
        WRITE(60, *) time
        WRITE(60, *) 3
        WRITE(60, *) t
        WRITE(60, *) 1
        WRITE(60, *) self%numNodes
        DO n=1, self%numNodes
          WRITE(60, "(I6,3(ES20.6E3))") n, output(n)%temperature
        END DO
        WRITE(60, "(A)") '$EndNodeData'
        CLOSE (60)

      END DO

    END SUBROUTINE printOutputCyl

    SUBROUTINE printCollisionsCyl(self, t)
      USE moduleRefParam
      USE moduleCaseParam
      USE moduleOutput
      IMPLICIT NONE

      CLASS(meshCylGeneric), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: n
      REAL(8):: time
      CHARACTER(:), ALLOCATABLE:: fileName
      CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations


      IF (collOutput) THEN
        time = DBLE(t)*tau*ti_ref
        WRITE(tstring, '(I6.6)') t

        fileName='OUTPUT_' // tstring// '_Collisions.msh'
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (60, file = path // folder // '/' //  fileName)
        WRITE(60, "(A)") '$MeshFormat'
        WRITE(60, "(A)") '2.2 0 8'
        WRITE(60, "(A)") '$EndMeshFormat'
        WRITE(60, "(A)") '$ElementData'
        WRITE(60, "(A)") '1'
        WRITE(60, "(A)") '"Collisions"'
        WRITE(60, *) 1
        WRITE(60, *) time
        WRITE(60, *) 3
        WRITE(60, *) t
        WRITE(60, *) 1
        WRITE(60, *) self%numVols
        DO n=1, self%numVols
          WRITE(60, "(I6,I10)") n + self%numEdges, self%vols(n)%obj%nColl
        END DO
        WRITE(60, "(A)") '$EndElementData'

        CLOSE(60)

      END IF

    END SUBROUTINE printCollisionsCyl

    SUBROUTINE printEMCyl(self, t)
      USE moduleRefParam
      USE moduleCaseParam
      USE moduleOutput
      IMPLICIT NONE

      CLASS(meshCylGeneric), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: n, e
      REAL(8):: time
      CHARACTER(:), ALLOCATABLE:: fileName
      CHARACTER (LEN=6):: tstring !TODO: Review to allow any number of iterations
      REAL(8):: xi(1:3)

      IF (emOutput) THEN
        time = DBLE(t)*tau*ti_ref
        WRITE(tstring, '(I6.6)') t

        fileName='OUTPUT_' // tstring// '_EMField.msh'
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (20, file = path // folder // '/' //  fileName)
        WRITE(20, "(A)") '$MeshFormat'
        WRITE(20, "(A)") '2.2 0 8'
        WRITE(20, "(A)") '$EndMeshFormat'
        WRITE(20, "(A)") '$NodeData'
        WRITE(20, "(A)") '1'
        WRITE(20, "(A)") '"Potential (V)"'
        WRITE(20, *) 1
        WRITE(20, *) time
        WRITE(20, *) 3
        WRITE(20, *) t
        WRITE(20, *) 1
        WRITE(20, *) self%numNodes
        DO n=1, self%numNodes
          WRITE(20, *) n, self%nodes(n)%obj%emData%phi*Volt_ref
        END DO
        WRITE(20, "(A)") '$EndNodeData'
        
        WRITE(20, "(A)") '$ElementData'
        WRITE(20, "(A)") '1'
        WRITE(20, "(A)") '"Electric Field (V/m)"'
        WRITE(20, *) 1
        WRITE(20, *) time
        WRITE(20, *) 3
        WRITE(20, *) t
        WRITE(20, *) 3
        WRITE(20, *) self%numVols
        DO e=1, self%numVols
          SELECT TYPE(elem=>self%vols(e)%obj)
          TYPE IS(meshVolCylQuad)
            xi = (/ 0.D0, 0.D0, 0.D0 /)
          TYPE IS(meshVolCylTria)
            xi = (/ 1.D0/3.D0, 1.D0/3.D0, 0.D0 /)
          END SELECT
          WRITE(20, *) e+self%numEdges, self%vols(e)%obj%gatherEF(xi)*EF_ref
        END DO
        WRITE(20, "(A)") '$EndElementData'
        CLOSE(20)

      END IF

    END SUBROUTINE printEMCyl
  
END MODULE moduleMeshCylRead