New functions to read VTU

Added a few functions to read .vtu meshes.

Output in .vtu changed to output data in 6 columns (seems to be the
    standard.)

src/modules/init/moduleInput.f90

@@ -860,6 +860,7 @@ MODULE moduleInput
     SUBROUTINE readGeometry(config)
       USE moduleMesh
       USE moduleMeshInputGmsh2, ONLY: initGmsh2
+      USE moduleMeshInputVTU,   ONLY: initVTU, readVTU !TEMPORARY TO TEST VTU OUTPUT
       USE moduleMeshInput0D,    ONLY: init0D
       USE moduleErrors
       USE moduleOutput
@@ -873,6 +874,7 @@ MODULE moduleInput
       LOGICAL:: found
       CHARACTER(:), ALLOCATABLE:: meshFormat, meshFile
       REAL(8):: volume
+      CHARACTER(:), ALLOCATABLE:: meshFileVTU !Temporary to test VTU OUTPUT
 
       object = 'geometry'
 
@@ -973,6 +975,8 @@ MODULE moduleInput
         CALL config%get(object // '.meshFile', meshFile, found)
         pathMeshParticle = path // meshFile
         CALL mesh%readMesh(pathMeshParticle)
+        meshFileVTU = '/home/jorge/fpakc/runs/cylFlow/mesh.vtu' !TEMPORARY TO TEST VTU OUTPUT
+        CALL readVTU(mesh, meshFileVTU) !TEMPORARY TO TEST VTU OUTPUT
         DEALLOCATE(meshFile)
         IF (doubleMesh) THEN
           !Reads the mesh file for collisions

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

@@ -1,8 +1,138 @@
-MODULE moduleMeshInputVTK
+MODULE moduleMeshInputVTU
   !Reads mesh in the VTU format
 
+  INTERFACE getValueFromLine
+    MODULE PROCEDURE getIntegerFromLine, getRealFromLine
+
+  END INTERFACE
+
   CONTAINS
-    SUBROUTINE initVTK(self)
+    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
+        PRINT *, iStart, iEnd
+        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
@@ -20,7 +150,7 @@ MODULE moduleMeshInputVTK
       END SELECT
       self%readMesh    => readVTU
 
-    END SUBROUTINE initVTK
+    END SUBROUTINE initVTU
 
     SUBROUTINE readVTU(self, filename)
       USE moduleMesh3DCart
@@ -34,6 +164,34 @@ MODULE moduleMeshInputVTK
       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
+
+      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 readIntegerBlock(fileID, numElements, entitiesID)
+      REWIND(fileID)
+
+      !Get the offsets to read connectivity
+      line = findLine(fileID, 'Name="offsets"')
+      ALLOCATE(offsets(1:numElements))
+      CALL readIntegerBlock(fileID, numElements, offsets)
+      REWIND(fileID)
+
+      CLOSE(fileID)
 
     END SUBROUTINE readVTU
 
@@ -47,4 +205,4 @@ MODULE moduleMeshInputVTK
 
     END SUBROUTINE readInitialVTU
 
-END MODULE moduleMeshInputVTK
+END MODULE moduleMeshInputVTU

src/modules/mesh/inout/vtu/moduleMeshOutputVTU.f90

@@ -103,45 +103,37 @@ MODULE moduleMeshOutputVTU
 
       CLASS(meshGeneric), INTENT(in):: self
       INTEGER, INTENT(in):: fileID
-      CHARACTER(LEN=25):: nodeFormat
-      CHARACTER(LEN=25):: cellFormat
       INTEGER:: e
-      INTEGER:: offset
+      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(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
+      WRITE(fileID, "(6(ES20.6E3))") (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes)
-      WRITE(fileID, nodeFormat) (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">'
-      DO e = 1, self%numCells
+      WRITE(fileID, "(6(I10))") (self%cells(e)%obj%getNodes(self%cells(e)%obj%nNodes) - 1, e = 1, self%numCells) !Array starts on 0
-        WRITE(cellFormat, "(A,I1,A)") "(",self%cells(e)%obj%nNodes,"(I10))"
-        WRITE(fileID, cellFormat, advance="no") self%cells(e)%obj%getNodes(self%cells(e)%obj%nNodes) - 1 !Array starts on 0
-
-      END DO
       WRITE(fileID, "(10X, A)") '</DataArray>'
       !Write offset of each cell
-      offset = 0
+      offset(1) = self%cells(1)%obj%nNodes
       WRITE(fileID, "(10X,A)") '<DataArray type="Int32" Name="offsets"  format="ascii">'
-      DO e = 1, self%numCells
+      DO e = 2, self%numCells
-        WRITE(cellFormat, "(A,I1,A)") "(I10)"
+        offset(e) = offset(e - 1) + self%cells(e)%obj%nNodes
-        offset = offset + self%cells(e)%obj%nNodes
-        WRITE(fileID, cellFormat, advance="no") offset
 
       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
-        WRITE(cellFormat, "(A,I1,A)") "(I10)"
+        types(e) = getCellType(self%cells(e)%obj)
-        WRITE(fileID, cellFormat, advance="no") getCellType(self%cells(e)%obj)
 
       END DO
+      WRITE(fileID, "(6(I10))") types
       WRITE(fileID, "(10X, A)") '</DataArray>'
       WRITE(fileID, "(8X, A)") '</Cells>'
 
@@ -156,7 +148,6 @@ MODULE moduleMeshOutputVTU
       INTEGER, INTENT(in):: fileID
       INTEGER, INTENT(in):: speciesIndex
       TYPE(outputFormat):: output(1:self%numNodes)
-      CHARACTER(LEN=25):: nodeFormat
       INTEGER:: n
 
       DO n = 1, self%numNodes
@@ -168,23 +159,19 @@ MODULE moduleMeshOutputVTU
       WRITE(fileID,"(8X,A)") '<PointData>'
       !Write density
       WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Density (m^-3)" NumberOfComponents="1">'
-      WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileID,"(6(ES20.6E3))") (output(n)%density, n  = 1, self%numNodes)
-      WRITE(fileID,nodeFormat) (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(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
+      WRITE(fileID,"(6(ES20.6E3))") (output(n)%velocity, n  = 1, self%numNodes)
-      WRITE(fileID,nodeFormat) (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(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileID,"(6(ES20.6E3))") (output(n)%pressure, n  = 1, self%numNodes)
-      WRITE(fileID,nodeFormat) (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(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileID,"(6(ES20.6E3))") (output(n)%temperature, n  = 1, self%numNodes)
-      WRITE(fileID,nodeFormat) (output(n)%temperature, n  = 1, self%numNodes)
       WRITE(fileID,"(10X,A)") '</DataArray>'
       !End node data
       WRITE(fileID,"(8X,A)") '</PointData>'
@@ -201,7 +188,6 @@ MODULE moduleMeshOutputVTU
       INTEGER:: k, c, n
       CHARACTER(:), ALLOCATABLE:: title
       CHARACTER (LEN=2):: cString
-      CHARACTER(LEN=25):: cellFormat
 
       WRITE(fileID,"(8X,A)") '<CellData>'
       DO k = 1, nCollPairs
@@ -209,8 +195,7 @@ MODULE moduleMeshOutputVTU
           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(cellFormat, "(A,I10, A)") "(", self%numCells, "(I10))"
+          WRITE(fileID, "(6(I10))") (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells)
-          WRITE(fileID, cellFormat) (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells)
           WRITE(fileID, "(10X, A)") '</DataArray>'
 
         END DO
@@ -226,33 +211,27 @@ MODULE moduleMeshOutputVTU
 
       CLASS(meshParticles), INTENT(in):: self
       INTEGER, INTENT(in):: fileID
-      CHARACTER(LEN=25):: nodeFormat
-      CHARACTER(LEN=25):: cellFormat
       INTEGER:: n
       REAL(8):: Xi(1:3)
 
-      Xi = (/ 0.D0, 0.D0, 0.D0 /)
-
       !Points in nodes
       WRITE(fileID,"(8X,A)") '<PointData>'
       !Electric potential
       WRITE(fileID,"(10X,A)") '<DataArray type="Float64" Name="Potential (V)" NumberOfComponents="1">'
-      WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%phi*Volt_ref, n  = 1, self%numNodes)
-      WRITE(fileID,nodeFormat) (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(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
+      WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%B*B_ref, n  = 1, self%numNodes)
-      WRITE(fileID,nodeFormat) (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(cellFormat, "(A,I10, A)") "(", self%numCells, "(3(ES20.6E3)))"
+      WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
-      WRITE(fileID, cellFormat) (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells)
       WRITE(fileID,"(10X,A)") '</DataArray>'
       WRITE(fileID,"(8X,A)") '</CellData>'
 
@@ -303,7 +282,6 @@ MODULE moduleMeshOutputVTU
       INTEGER, INTENT(in):: speciesIndex
       TYPE(outputFormat):: outputMean(1:self%numNodes)
       TYPE(outputFormat):: outputDeviation(1:self%numNodes)
-      CHARACTER(LEN=25):: nodeFormat
       INTEGER:: n
 
       DO n = 1, self%numNodes
@@ -320,33 +298,29 @@ MODULE moduleMeshOutputVTU
       !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(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileIDMean,     "(6(ES20.6E3))") (outputMean(n)%density,      n  = 1, self%numNodes)
-      WRITE(fileIDMean,     nodeFormat) (outputMean(n)%density,      n  = 1, self%numNodes)
+      WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%density, n  = 1, self%numNodes)
-      WRITE(fileIDDeviation,nodeFormat) (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(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))"
+      WRITE(fileIDMean,     "(6(ES20.6E3))") (outputMean(n)%velocity,      n  = 1, self%numNodes)
-      WRITE(fileIDMean,     nodeFormat) (outputMean(n)%velocity,      n  = 1, self%numNodes)
+      WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%velocity, n  = 1, self%numNodes)
-      WRITE(fileIDDeviation,nodeFormat) (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(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileIDMean,     "(6(ES20.6E3))") (outputMean(n)%pressure,      n  = 1, self%numNodes)
-      WRITE(fileIDMean,     nodeFormat) (outputMean(n)%pressure,      n  = 1, self%numNodes)
+      WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%pressure, n  = 1, self%numNodes)
-      WRITE(fileIDDeviation,nodeFormat) (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(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))"
+      WRITE(fileIDMean,     "(6(ES20.6E3))") (outputMean(n)%temperature,      n  = 1, self%numNodes)
-      WRITE(fileIDMean,     nodeFormat) (outputMean(n)%temperature,      n  = 1, self%numNodes)
+      WRITE(fileIDDeviation,"(6(ES20.6E3))") (outputDeviation(n)%temperature, n  = 1, self%numNodes)
-      WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%temperature, n  = 1, self%numNodes)
       WRITE(fileIDMean,     "(10X,A)") '</DataArray>'
       WRITE(fileIDDeviation,"(10X,A)") '</DataArray>'
       !End node data