MODULE moduleMeshCylRead
  USE moduleMesh
  USE moduleMeshCyl
  USE moduleMeshCylBoundary

  TYPE, EXTENDS(meshGeneric):: meshCylGeneric
    CONTAINS
      PROCEDURE, PASS:: readMesh    => readMeshCyl

  END TYPE

  INTERFACE connected
    MODULE PROCEDURE connectedVolVol, connectedVolEdge

  END INTERFACE connected

  CONTAINS
    SUBROUTINE readMeshCyl(self, filename)
      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.
        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
        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

    !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

END MODULE moduleMeshCylRead