MODULE moduleMeshInputVTU
  !Reads mesh in the VTU format

  INTERFACE getValueFromLine
    MODULE PROCEDURE getIntegerFromLine, getRealFromLine

  END INTERFACE

  INTERFACE readDataBlock
    MODULE PROCEDURE readIntegerBlock, readRealBlock

  END INTERFACE

  CONTAINS
    FUNCTION findLine(fileID, text) RESULT(line)
      USE moduleErrors
      IMPLICIT NONE
      
      INTEGER, INTENT(in):: fileID
      CHARACTER(*):: text
      CHARACTER(LEN=256):: line
      INTEGER:: error, found

      error = 0
      found = 0
      !Reads all the file to find the 'text' string in a line
      DO WHILE(error == 0 .AND. found == 0)
        READ(fileID, "(A)", IOSTAT=error) line
        found = INDEX(line, text)
        IF (found > 0) THEN
          EXIT

        END IF
      END DO

      !If no line is found, return an error
      IF (found == 0) THEN
        CALL criticalError('String ' // text // ' not found in file.', 'findLine')

      END IF

    END FUNCTION findLine

    SUBROUTINE getIntegerFromLine(line, label, valueInteger)
      USE moduleErrors
      IMPLICIT NONE

      CHARACTER(LEN=256), INTENT(in):: line
      CHARACTER(*), INTENT(in):: label
      INTEGER, INTENT(out):: valueInteger
      INTEGER:: labelStart, valueStart, valueEnd

      labelStart = 0
      labelStart = INDEX(line, label)
      IF (labelStart == 0) THEN
        CALL criticalError('Label '// label // ' not found in line', 'getValueFromLine')

      END IF
      valueStart = INDEX(line(labelStart:256), '"') + labelStart
      valueEnd   = INDEX(line(valueStart:256), '"') + valueStart - 2
      READ(line(valueStart:valueEnd), '(I10)') valueInteger

    END SUBROUTINE getIntegerFromLine
    
    SUBROUTINE getRealFromLine(line, label, valueReal)
      USE moduleErrors
      IMPLICIT NONE

      CHARACTER(LEN=256), INTENT(in):: line
      CHARACTER(*), INTENT(in):: label
      REAL(8), INTENT(out):: valueReal
      INTEGER:: labelStart, valueStart, valueEnd

      labelStart = 0
      labelStart = INDEX(line, label)
      IF (labelStart == 0) THEN
        CALL criticalError('Label '// label // ' not found in line', 'getValueFromLine')

      END IF
      valueStart = INDEX(line(labelStart:256), '"') + labelStart
      valueEnd   = INDEX(line(valueStart:256), '"') + valueStart - 2
      READ(line(valueStart:valueEnd), '(F16.14)') valueReal

    END SUBROUTINE getRealFromLine

    SUBROUTINE readIntegerBlock(fileID, nData, array)
      IMPLICIT NONE

      INTEGER, INTENT(in):: fileID
      INTEGER, INTENT(in):: nData
      INTEGER, INTENT(out):: array(1:nData)
      INTEGER:: iStart, iEnd, block

      iStart = 0
      iEnd = 0
      block = 6 !Assumes block of data in 6 columns

      DO WHILE (iStart < nData)
        iStart = iStart + 1
        iEnd   = iStart - 1 + block
        IF (iEnd > nData) THEN
          iEnd = nData

        END IF
        READ(fileID, *) array(iStart:iEnd)
        iStart = iEnd

      END DO

    END SUBROUTINE readIntegerBlock

    SUBROUTINE readRealBlock(fileID, nData, array)
      IMPLICIT NONE

      INTEGER, INTENT(in):: fileID
      INTEGER, INTENT(in):: nData
      REAL(8), INTENT(out):: array(1:nData)
      INTEGER:: iStart, iEnd, block

      iStart = 0
      iEnd = 0
      block = 6 !Assumes block of data in 6 columns

      DO WHILE (iStart < nData)
        iStart = iStart + 1
        iEnd   = iStart - 1 + block
        IF (iEnd > nData) THEN
          iEnd = nData

        END IF
        READ(fileID, *) array(iStart:iEnd)
        iStart = iEnd

      END DO

    END SUBROUTINE readRealBlock

    SUBROUTINE initVTU(self)
      USE moduleMesh
      USE moduleMeshOutputVTU
      IMPLICIT NONE

      CLASS(meshGeneric), INTENT(inout), TARGET:: self

      IF (ASSOCIATED(meshForMCC, self)) self%printColl => printCollVTU
      SELECT TYPE(self)
      TYPE IS(meshParticles)
        self%printOutput  => printOutputVTU
        self%printEM      => printEMVTU
        self%readInitial  => readInitialVTU
        self%printAverage => printAverageVTU

      END SELECT
      self%readMesh    => readVTU

    END SUBROUTINE initVTU

    SUBROUTINE readVTU(self, filename)
      USE moduleMesh3DCart
      USE moduleMesh2DCyl
      USE moduleMesh2DCart
      USE moduleMesh1DRad
      USE moduleMesh1DCart
      IMPLICIT NONE

      CLASS(meshGeneric), INTENT(inout):: self
      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
      REAL(8):: r(1:3) !3 generic coordinates
      INTEGER:: fileID
      CHARACTER(LEN=256):: line
      INTEGER:: numNodes, numElements, numEdges
      INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types
      REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates
      INTEGER:: n, e, c
      INTEGER, ALLOCATABLE:: p(:)

      fileID = 10

      OPEN(fileID, FILE=TRIM(filename))

      !Find the number of nodes and elements (edges+cells) in the mesh.
      line = findLine(fileID, '<Piece')
      CALL getValueFromLine(line, 'NumberOfPoints', numNodes)
      CALL getValueFromLine(line, 'NumberOfCells', numElements)
      REWIND(fileID)

      !Get the IDs of the cells to identify physical surfaces
      line = findLine(fileID, 'Name="CellEntityIds"')
      ALLOCATE(entitiesID(1:numElements))
      CALL readDataBlock(fileID, numElements, entitiesID)
      REWIND(fileID)

      !Get the offsets to read connectivity
      line = findLine(fileID, 'Name="offsets"')
      ALLOCATE(offsets(1:numElements))
      CALL readDataBlock(fileID, numElements, offsets)
      REWIND(fileID)

      !Get the connectivity of elements to nodes
      line = findline(fileID, 'Name="connectivity"')
      ALLOCATE(connectivity(1:MAXVAL(offsets)))
      CALL readDataBlock(fileID, MAXVAL(offsets), connectivity)
      !Shift connectivity to start in 1
      connectivity = connectivity + 1
      REWIND(fileID)

      !Get the type of elements
      line = findline(fileID, 'Name="types"')
      ALLOCATE(types(1:numElements))
      CALL readDataBlock(fileID, numElements, types)
      REWIND(fileID)

      !Get nodes coordinates
      line = findline(fileID, 'Name="Points"')
      ALLOCATE(coordinates(1:3*numNodes))
      CALL readDataBlock(fileID, 3*numNodes, coordinates)
      REWIND(fileID)

      CLOSE(fileID)

      !All relevant information from the .vtu file has been read. Time to build the mesh.
      self%numNodes = numNodes
      ALLOCATE(self%nodes(1:self%numNodes))

      SELECT TYPE(self)
      TYPE IS(meshParticles)
        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

      END SELECT

      DO n = 1, self%numNodes
        !Get the coordinates for each direction
        r(1) = coordinates(3*(n-1)+1)
        r(2) = coordinates(3*(n-1)+2)
        r(3) = coordinates(3*(n-1)+3)

        !Allocate node
        SELECT CASE(self%dimen)
        CASE(3)
          ALLOCATE(meshNode3Dcart::self%nodes(n)%obj)

        CASE(2)
          SELECT CASE(self%geometry)
          CASE("Cyl")
            ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj)

          CASE("Cart")
            ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)

          END SELECT

          r(3) = 0.D0

        CASE(1)
          SELECT CASE(self%geometry)
          CASE("Rad")
            ALLOCATE(meshNode1DRad:: self%nodes(n)%obj)

          CASE("Cart")
            ALLOCATE(meshNode1DCart:: self%nodes(n)%obj)

          END SELECT
          r(2:3) = 0.D0

        END SELECT

        !Init node
        CALL self%nodes(n)%obj%init(n, r)

      END DO

      !Count the number of edges
      numEdges = 0
      SELECT CASE(self%dimen)
      CASE(3)
        !Edges are triangles, type 5 in VTK
        numEdges = COUNT(types==5)
      
      CASE(2)
        !Edges are segments, type 3 in VTK
        numEdges = COUNT(types==3)

      CASE(1)
        !Edges are nodes, type 1 in VTK
        numEdges = COUNT(types==1)
      
      END SELECT

      self%numCells = numElements - numEdges

      SELECT TYPE(self)
      TYPE IS(meshParticles)

        self%numEdges = numEdges

        !Allocate array of edges
        ALLOCATE(self%edges(1:self%numEdges))

      END SELECT

      !Allocates array of cells
      ALLOCATE(self%cells(1:self%numCells))

      !Read edges only if mesh is for tracking particles, not for collisions
      e = 0
      SELECT TYPE(self)
      TYPE IS(meshParticles)
        DO n = 1, numElements
          SELECT CASE(self%dimen)
          CASE(3)
            IF (types(n) == 5) THEN
              e = e + 1

              ALLOCATE(p(1:3))
              p(1) = connectivity(offsets(n) - 2)
              p(2) = connectivity(offsets(n) - 1)
              p(3) = connectivity(offsets(n))

              !Allocate edge
              ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj)

            END IF

          CASE(2)
            IF (types(n) == 3) THEN
              e = e + 1
              ALLOCATE(p(1:2))
              p(1) = connectivity(offsets(n) - 1)
              p(2) = connectivity(offsets(n))

              !Allocate edge
              SELECT CASE(self%geometry)
              CASE("Cyl")
                ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)

              CASE("Cart")
                ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)

              END SELECT

            END IF

          CASE(1)
            IF (types(n) == 1) THEN
              e = e + 1
              ALLOCATE(p(1:1))
              p(1) = connectivity(offsets(n))

              !Allocate edge
              SELECT CASE(self%geometry)
              CASE("Rad")
                ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)

              CASE("Cart")
                ALLOCATE(meshEdge1DCart:: self%edges(e)%obj)

              END SELECT

            END IF

          END SELECT

          ! If an edge was found, init it
          if (allocated(p)) then
            call self%edges(e)%obj%init(n, p, entitiesID(n))

            deallocate(p)

          end if

        END DO

      END SELECT

      !Read cells
      c = 0
      DO n = 1, numElements
        SELECT CASE(self%dimen)
        CASE(3)
          SELECT CASE(types(n))
          CASE(10)
            !Thetraedron
            c = c + 1

            ALLOCATE(p(1:4))
            p(1) = connectivity(offsets(n) - 3)
            p(2) = connectivity(offsets(n) - 2)
            p(3) = connectivity(offsets(n) - 1)
            p(4) = connectivity(offsets(n))

            !Allocate cell
            ALLOCATE(meshCell3DCartTetra::  self%cells(c)%obj)

            !Init cell
            CALL self%cells(c)%obj%init(c, p, self%nodes)
            DEALLOCATE(p)

          END SELECT

        CASE(2)
          SELECT CASE(types(n))
          CASE(5)
            !Triangular element
            c = c + 1

            ALLOCATE(p(1:3))
            p(1) = connectivity(offsets(n) - 2)
            p(2) = connectivity(offsets(n) - 1)
            p(3) = connectivity(offsets(n))

            !Allocate cell
            SELECT CASE(self%geometry)
            CASE("Cyl")
              ALLOCATE(meshCell2DCylTria:: self%cells(c)%obj)

            CASE("Cart")
              ALLOCATE(meshCell2DCartTria:: self%cells(c)%obj)

            END SELECT

            !Init cell
            CALL self%cells(c)%obj%init(c, p, self%nodes)
            DEALLOCATE(p)

          CASE(9)
            !Quadrilateral element
            c = c + 1

            ALLOCATE(p(1:4))
            p(1) = connectivity(offsets(n) - 3)
            p(2) = connectivity(offsets(n) - 2)
            p(3) = connectivity(offsets(n) - 1)
            p(4) = connectivity(offsets(n))

            !Allocate cell
            SELECT CASE(self%geometry)
            CASE("Cyl")
              ALLOCATE(meshCell2DCylQuad:: self%cells(c)%obj)

            CASE("Cart")
              ALLOCATE(meshCell2DCartQuad:: self%cells(c)%obj)

            END SELECT

            !Init cell
            CALL self%cells(c)%obj%init(c, p, self%nodes)
            DEALLOCATE(p)


          END SELECT

        CASE(1)
          SELECT CASE(types(n))
          CASE(3)
            !Segment element
            c = c + 1

            ALLOCATE(p(1:2))
            p(1) = connectivity(offsets(n) - 1)
            p(2) = connectivity(offsets(n))

            !Allocate cell
            SELECT CASE(self%geometry)
            CASE("Rad")
              ALLOCATE(meshCell1DRadSegm:: self%cells(c)%obj)

            CASE("Cart")
              ALLOCATE(meshCell1DCartSegm:: self%cells(c)%obj)

            END SELECT

            !Init cell
            CALL self%cells(c)%obj%init(c, p, self%nodes)
            DEALLOCATE(p)

          END SELECT

        END SELECT

      END DO

      !Call mesh connectivity
      CALL self%connectMesh

    END SUBROUTINE readVTU

    SUBROUTINE readInitialVTU(filename, density, velocity, temperature)
      IMPLICIT NONE

      CHARACTER(:), ALLOCATABLE, INTENT(in):: filename
      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: density
      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:,:):: velocity
      REAL(8), ALLOCATABLE, INTENT(out), DIMENSION(:):: temperature
      INTEGER:: fileID
      CHARACTER(:), ALLOCATABLE:: line
      INTEGER:: numNodes
      REAL(8), ALLOCATABLE:: velocityBlock(:)
      INTEGER:: n

      fileID = 10

      OPEN(fileID, file = TRIM(filename))

      line = findLine(fileID, '<Piece')
      CALL getValueFromLine(line, 'NumberOfPoints', numNodes)

      !Read the species density
      line = findLine(fileID, 'Name="Density')
      ALLOCATE(density(1:numNodes))
      CALL readDataBlock(fileID, numNodes, density)
      REWIND(fileID)

      !Read the species velocity
      line = findLine(fileID, 'Name="Velocity')
      ALLOCATE(velocityBlock(1:3*numNodes))
      CALL readDataBlock(fileID, 3*numNodes, velocityBlock)
      ALLOCATE(velocity(1:numNodes, 1:3))
      DO n = 1, numNodes
        velocity(n, 1) = velocityBlock(3*(n-1)+1)
        velocity(n, 2) = velocityBlock(3*(n-1)+2)
        velocity(n, 3) = velocityBlock(3*(n-1)+3)

      END DO
      DEALLOCATE(velocityBlock)
      REWIND(fileID)

      !Read the species temperature
      line = findLine(fileID, 'Name="Temperature')
      ALLOCATE(temperature(1:numNodes))
      CALL readDataBlock(fileID, numNodes, temperature)
      REWIND(fileID)

      close(fileID)

    END SUBROUTINE readInitialVTU

END MODULE moduleMeshInputVTU