fpakc/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90

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
        PRINT *, iStart, iEnd
        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
      USE moduleBoundary
      IMPLICIT NONE

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

      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)
      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.

    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

    END SUBROUTINE readInitialVTU

END MODULE moduleMeshInputVTU