MODULE moduleMeshOutputVTU

  CONTAINS

    SUBROUTINE writeHeader(nNodes, nCells, fileID)
      USE moduleMesh
      IMPLICIT NONE

      INTEGER, INTENT(in):: nNodes, nCells
      INTEGER, INTENT(in):: fileID

      WRITE(fileID,"(A)") '<?xml version="1.0"?>'
      WRITE(fileID,"(2X, A)") '<VTKFile type="UnstructuredGrid">'
      WRITE(fileID,"(4X, A,ES20.6E3,A)") '<UnstructuredGrid>'
      WRITE(fileID,"(6X, A, I10, A, I10, A)") '<Piece NumberOfPoints="', nNodes, '" NumberOfCells="', nCells, '">'

    END SUBROUTINE writeHeader

    SUBROUTINE writeFooter(fileID)
      IMPLICIT NONE

      INTEGER, INTENT(in):: fileID

      WRITE(fileID,"(6X, A)") '</Piece>'
      WRITE(fileID,"(4X, A)") '</UnstructuredGrid>'
      WRITE(fileID,"(2X, A)") '</VTKFile>'

    END SUBROUTINE writeFooter

    FUNCTION getCellType(cell) RESULT(indexType)
      USE moduleMesh3DCart
      USE moduleMesh2DCyl
      USE moduleMesh2DCart
      USE moduleMesh1DRad
      USE moduleMesh1DCart
      USE moduleMesh0D
      USE moduleErrors
      IMPLICIT NONE

      CLASS(meshCell), INTENT(in):: cell
      INTEGER:: indexType
      
      indexType = 0

      SELECT TYPE(cell)
      TYPE IS(meshCell3DCartTetra)
        indexType = 10

      TYPE IS(meshCell2DCylQuad)
        indexType = 9

      TYPE IS(meshCell2DCartQuad)
        indexType = 9

      TYPE IS(meshCell2DCylTria)
        indexType = 5

      TYPE IS(meshCell2DCartTria)
        indexType = 5

      TYPE IS(meshCell1DRadSegm)
        indexType = 3

      TYPE IS(meshCell1DCartSegm)
        indexType = 3

      TYPE IS(meshCell0D)
        indexType = 1

      CLASS DEFAULT
        CALL criticalError('Cell not valid for VTU output', 'getCellType')

      END SELECT

    END FUNCTION getCellType

    SUBROUTINE writeMesh(self, fileID)
      USE moduleMesh
      USE moduleRefParam
      IMPLICIT NONE

      CLASS(meshGeneric), INTENT(in):: self
      INTEGER, INTENT(in):: fileID
      INTEGER:: e
      INTEGER:: offset(1:self%numCells), types(1:self%numCells)

      !Write nodes coordinates
      WRITE(fileID, "(8X, A)") '<Points>'
      WRITE(fileID, "(10X,A)") '<DataArray type="Float32" NumberOfComponents="3" format="ascii">'
      WRITE(fileID, "(6(ES20.6E3))") (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes)
      WRITE(fileID, "(10X, A)") '</DataArray>'
      WRITE(fileID, "(8X, A)") '</Points>'

      WRITE(fileID, "(8X, A)") '<Cells>'
      !Write nodes connectivity of each cell
      WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="connectivity"  format="ascii">'
      WRITE(fileID, "(6(I10))") (self%cells(e)%obj%getNodes(self%cells(e)%obj%nNodes) - 1, e = 1, self%numCells) !Array starts on 0
      WRITE(fileID, "(10X, A)") '</DataArray>'
      !Write offset of each cell
      offset(1) = self%cells(1)%obj%nNodes
      WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="offsets"  format="ascii">'
      DO e = 2, self%numCells
        offset(e) = offset(e - 1) + self%cells(e)%obj%nNodes

      END DO
      WRITE(fileID, "(6(I10))") offset
      WRITE(fileID, "(10X, A)") '</DataArray>'
      !Write type of each cell
      WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="types"  format="ascii">'
      DO e = 1, self%numCells
        types(e) = getCellType(self%cells(e)%obj)

      END DO
      WRITE(fileID, "(6(I10))") types
      WRITE(fileID, "(10X, A)") '</DataArray>'
      WRITE(fileID, "(8X, A)") '</Cells>'

    END SUBROUTINE writeMesh

    SUBROUTINE writeSpeciesOutput(self, fileID, speciesIndex)
      USE moduleMesh
      USE moduleOutput
      IMPLICIT NONE

      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: fileID
      INTEGER, INTENT(in):: speciesIndex
      TYPE(outputFormat):: output(1:self%numNodes)
      INTEGER:: n

      DO n = 1, self%numNodes
        CALL calculateOutput(self%nodes(n)%obj%output(speciesIndex), output(n), self%nodes(n)%obj%v, species(speciesIndex)%obj)

      END DO

      !Write node data
      WRITE(fileID,"(8X,A)") '<PointData>'
      !Write density
      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Density (m^-3)" NumberOfComponents="1">'
      WRITE(fileID,"(6(ES20.6E3))") (output(n)%density, n  = 1, self%numNodes)
      WRITE(fileID,"(10X,A)") '</DataArray>'
      !Write velocity
      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Velocity (m s^-1)" NumberOfComponents="3">'
      WRITE(fileID,"(6(ES20.6E3))") (output(n)%velocity, n  = 1, self%numNodes)
      WRITE(fileID,"(10X,A)") '</DataArray>'
      !Write pressure
      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Pressure (Pa)" NumberOfComponents="1">'
      WRITE(fileID,"(6(ES20.6E3))") (output(n)%pressure, n  = 1, self%numNodes)
      WRITE(fileID,"(10X,A)") '</DataArray>'
      !Write temperature
      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Temperature (K)" NumberOfComponents="1">'
      WRITE(fileID,"(6(ES20.6E3))") (output(n)%temperature, n  = 1, self%numNodes)
      WRITE(fileID,"(10X,A)") '</DataArray>'
      !End node data
      WRITE(fileID,"(8X,A)") '</PointData>'

    END SUBROUTINE writeSpeciesOutput

    SUBROUTINE writeCollOutput(self,fileID)
      USE moduleMesh
      USE moduleCollisions
      IMPLICIT NONE

      CLASS(meshGeneric), INTENT(in):: self
      INTEGER, INTENT(in):: fileID
      INTEGER:: k, c, n
      CHARACTER(:), ALLOCATABLE:: title
      CHARACTER (LEN=2):: cString

      WRITE(fileID,"(8X,A)") '<CellData>'
      DO k = 1, nCollPairs
        DO c = 1, interactionMatrix(k)%amount
          WRITE(cString, "(I2)") c
          title = 'Pair ' // interactionMatrix(k)%sp_i%name // '-' // interactionMatrix(k)%sp_j%name // ' collision ' // cString
          WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="',title, '" NumberOfComponents="1">'
          WRITE(fileID, "(6(I10))") (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells)
          WRITE(fileID, "(10X, A)") '</DataArray>'

        END DO
      END DO
      WRITE(fileID,"(8X,A)") '</CellData>'

    END SUBROUTINE writeCollOutput

    SUBROUTINE writeEM(self, fileID)
      USE moduleMesh
      USE moduleRefParam
      IMPLICIT NONE

      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: fileID
      INTEGER:: n
      REAL(8):: Xi(1:3)

      !Points in nodes
      WRITE(fileID,"(8X,A)") '<PointData>'
      !Electric potential
      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Potential (V)" NumberOfComponents="1">'
      WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%phi*Volt_ref, n  = 1, self%numNodes)
      WRITE(fileID,"(10X,A)") '</DataArray>'
      !Magnetic Field
      WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Magnetic Field (T)" NumberOfComponents="3">'
      WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%B*B_ref, n  = 1, self%numNodes)
      WRITE(fileID,"(10X,A)") '</DataArray>'
      WRITE(fileID,"(8X,A)") '</PointData>'

      !Cell Data
      Xi = (/ 0.D0, 0.D0, 0.D0 /)
      WRITE(fileID,"(8X,A)") '<CellData>'
      !Electric field
      WRITE(fileID,"(10X,A, A, A)") '<DataArray type="Float64" Name="Electric Field (V m^-1)" NumberOfComponents="3">'
      WRITE(fileID,"(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
      WRITE(fileID,"(10X,A)") '</DataArray>'
      WRITE(fileID,"(8X,A)") '</CellData>'

    END SUBROUTINE writeEM

    SUBROUTINE writeCollection(fileID, t, fileNameStep, fileNameCollection)
      USE moduleCaseParam
      USE moduleOutput
      USE moduleRefParam
      IMPLICIT NONE

      INTEGER:: fileID
      INTEGER, INTENT(in):: t
      CHARACTER(*):: fileNameStep, fileNameCollection

      IF (t == tInitial) THEN
        !Create collection file
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameCollection
        OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection)
        WRITE (fileID + 1, "(A)") '<VTKFile type="Collection">'
        WRITE (fileID + 1, "(2X, A)") '<Collection>'
        CLOSE(fileID + 1)

      END IF

      !Write iteration file in collection
      OPEN (fileID + 1, file = path // folder  // '/' // fileNameCollection, ACCESS='APPEND')
      WRITE(fileID + 1, "(4X, A, ES20.6E3, A, A, A)") '<DataSet timestep="', DBLE(t)*tauMin*ti_ref,'" file="', fileNameStep,'"/>'

      !Close collection file
      IF (t == tFinal) THEN
        WRITE (fileID + 1, "(2X, A)") '</Collection>'
        WRITE (fileID + 1, "(A)") '</VTKFile>'

      END IF
      CLOSE(fileID + 1)

    END SUBROUTINE writeCollection
    
    SUBROUTINE writeAverage(self, fileIDMean, fileIDDeviation, speciesIndex)
      USE moduleMesh
      USE moduleOutput
      USE moduleAverage
      IMPLICIT NONE

      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: fileIDMean, fileIDDeviation
      INTEGER, INTENT(in):: speciesIndex
      TYPE(outputFormat):: outputMean(1:self%numNodes)
      TYPE(outputFormat):: outputDeviation(1:self%numNodes)
      INTEGER:: n

      DO n = 1, self%numNodes
        CALL calculateOutput(averageScheme(n)%mean%output(speciesIndex),      outputMean(n),      &
                             self%nodes(n)%obj%v, species(speciesIndex)%obj)
        CALL calculateOutput(averageScheme(n)%deviation%output(speciesIndex), outputDeviation(n), &
                             self%nodes(n)%obj%v, species(speciesIndex)%obj)

      END DO

      !Write node data
      WRITE(fileIDMean,     "(8X,A)") '<PointData>'
      WRITE(fileIDDeviation,"(8X,A)") '<PointData>'
      !Write density
      WRITE(fileIDMean,     "(10X,A)") '<DataArray type="Float64" Name="Density, mean (m^-3)"      NumberOfComponents="1">'
      WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Density, deviation (m^-3)" NumberOfComponents="1">'
      WRITE(fileIDMean,     "(6(ES20.6E3))") (outputMean(n)%density,      n  = 1, self%numNodes)
      WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%density, n  = 1, self%numNodes)
      WRITE(fileIDMean,"(10X,A)")      '</DataArray>'
      WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
      !Write velocity
      WRITE(fileIDMean,     "(10X,A)") '<DataArray type="Float64" Name="Velocity, mean (m s^-1)"      NumberOfComponents="3">'
      WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Velocity, deviation (m s^-1)" NumberOfComponents="3">'
      WRITE(fileIDMean,     "(6(ES20.6E3))") (outputMean(n)%velocity,      n  = 1, self%numNodes)
      WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%velocity, n  = 1, self%numNodes)
      WRITE(fileIDMean,     "(10X,A)") '</DataArray>'
      WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
      !Write pressure
      WRITE(fileIDMean,     "(10X,A)") '<DataArray type="Float64" Name="Pressure, mean (Pa)"      NumberOfComponents="1">'
      WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Pressure, deviation (Pa)" NumberOfComponents="1">'
      WRITE(fileIDMean,     "(6(ES20.6E3))") (outputMean(n)%pressure,      n  = 1, self%numNodes)
      WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%pressure, n  = 1, self%numNodes)
      WRITE(fileIDMean,     "(10X,A)") '</DataArray>'
      WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
      !Write temperature
      WRITE(fileIDMean,     "(10X,A)") '<DataArray type="Float64" Name="Temperature, mean (K)"      NumberOfComponents="1">'
      WRITE(fileIDDeviation,"(10X,A)") '<DataArray type="Float64" Name="Temperature, deviation (K)" NumberOfComponents="1">'
      WRITE(fileIDMean,     "(6(ES20.6E3))") (outputMean(n)%temperature,      n  = 1, self%numNodes)
      WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%temperature, n  = 1, self%numNodes)
      WRITE(fileIDMean,     "(10X,A)") '</DataArray>'
      WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
      !End node data
      WRITE(fileIDMean,     "(8X,A)") '</PointData>'
      WRITE(fileIDDeviation,"(8X,A)") '</PointData>'

    END SUBROUTINE writeAverage

    SUBROUTINE printOutputVTU(self,t)
      USE moduleMesh
      USE moduleSpecies
      USE moduleMeshInoutCommon
      IMPLICIT NONE

      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: i, fileID
      CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection

      fileID = 60

      DO i = 1, nSpecies
        fileName = formatFileName(prefix, species(i)%obj%name, 'vtu', t)
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (fileID, file = path // folder // '/' //  fileName)

        CALL writeHeader(self%numNodes, self%numCells, fileID)

        CALL writeMesh(self, fileID)

        CALL writeSpeciesOutput(self, fileID, i)

        CALL writeFooter(fileID)

        CLOSE(fileID)

        !Write collection file for time plotting
        fileNameCollection = formatFileName('Collection', species(i)%obj%name, 'pvd')
        CALL writeCollection(fileID, t, fileName, filenameCollection)

      END DO

    END SUBROUTINE printOutputVTU

    SUBROUTINE printCollVTU(self,t)
      USE moduleMesh
      USE moduleOutput
      USE moduleMeshInoutCommon
      IMPLICIT NONE

      CLASS(meshGeneric), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: fileID
      CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection
      CHARACTER (LEN=iterationDigits):: tstring

      fileID = 62

      IF (collOutput) THEN
        fileName = formatFileName(prefix, 'Collisions', 'vtu', t)
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (fileID, file = path // folder // '/' //  fileName)

        CALL writeHeader(self%numNodes, self%numCells, fileID)

        CALL writeMesh(self, fileID)

        CALL writeCollOutput(self, fileID)

        CALL writeFooter(fileID)

        CLOSE(fileID)

        !Write collection file for time plotting
        fileNameCollection = formatFileName('Collection', 'Collisions', 'pvd')
        CALL writeCollection(fileID, t, fileName, filenameCollection)

      END IF

    END SUBROUTINE printCollVTU

    SUBROUTINE printEMVTU(self, t)
      USE moduleMesh
      USE moduleMeshInoutCommon
      IMPLICIT NONE

      CLASS(meshParticles), INTENT(in):: self
      INTEGER, INTENT(in):: t
      INTEGER:: fileID
      CHARACTER(:), ALLOCATABLE:: fileName, fileNameCollection

      fileID = 64

      IF (emOutput) THEN
        fileName = formatFileName(prefix, 'EMField', 'vtu', t)
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileName
        OPEN (fileID, file = path // folder // '/' //  fileName)

        CALL writeHeader(self%numNodes, self%numCells, fileID)

        CALL writeMesh(self, fileID)

        CALL writeEM(self, fileID)

        CALL writeFooter(fileID)

        CLOSE(fileID)

        !Write collection file for time plotting
        fileNameCollection = formatFileName('Collection', 'EMField', 'pvd')
        CALL writeCollection(fileID, t, fileName, filenameCollection)

      END IF

    END SUBROUTINE printEMVTU

    SUBROUTINE printAverageVTU(self)
      USE moduleMesh
      USE moduleSpecies
      USE moduleMeshInoutCommon
      IMPLICIT NONE

      CLASS(meshParticles), INTENT(in):: self
      INTEGER:: i, fileIDMean, fileIDDeviation
      CHARACTER(:), ALLOCATABLE:: fileNameMean, fileNameDeviation

      fileIDMean = 66
      fileIDDeviation = 67

      DO i = 1, nSpecies
        fileNameMean = formatFileName('Average_mean', species(i)%obj%name, 'vtu')
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameMean
        OPEN (fileIDMean, file = path // folder // '/' //  fileNameMean)

        fileNameDeviation = formatFileName('Average_deviation', species(i)%obj%name, 'vtu')
        WRITE(*, "(6X,A15,A)") "Creating file: ", fileNameDeviation
        OPEN (fileIDDeviation, file = path // folder // '/' //  fileNameDeviation)

        CALL writeHeader(self%numNodes, self%numCells, fileIDMean)
        CALL writeHeader(self%numNodes, self%numCells, fileIDDeviation)

        CALL writeMesh(self, fileIDMean)
        CALL writeMesh(self, fileIDDeviation)

        CALL writeAverage(self, fileIDMean, fileIDDeviation, i)

        CALL writeFooter(fileIDMean)
        CALL writeFooter(fileIDDeviation)

        CLOSE(fileIDMean)
        CLOSE(fileIDDeviation)

      END DO

    END SUBROUTINE printAverageVTU

END MODULE moduleMeshOutputVTU