diff --git a/src/modules/init/moduleInput.f90 b/src/modules/init/moduleInput.f90 index 0a33903..82d458e 100644 --- a/src/modules/init/moduleInput.f90 +++ b/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 diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 index da6c792..505f10d 100644 --- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 +++ b/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, '' WRITE(fileID, "(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" - WRITE(fileID, nodeFormat) (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes) + WRITE(fileID, "(6(ES20.6E3))") (self%nodes(e)%obj%getCoordinates()*L_ref, e = 1, self%numNodes) WRITE(fileID, "(10X, A)") '' WRITE(fileID, "(8X, A)") '' WRITE(fileID, "(8X, A)") '' !Write nodes connectivity of each cell WRITE(fileID, "(10X,A)") '' - DO e = 1, self%numCells - 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, "(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)") '' !Write offset of each cell - offset = 0 + offset(1) = self%cells(1)%obj%nNodes WRITE(fileID, "(10X,A)") '' - DO e = 1, self%numCells - WRITE(cellFormat, "(A,I1,A)") "(I10)" - offset = offset + self%cells(e)%obj%nNodes - WRITE(fileID, cellFormat, advance="no") offset + 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)") '' !Write type of each cell WRITE(fileID, "(10X,A)") '' DO e = 1, self%numCells - WRITE(cellFormat, "(A,I1,A)") "(I10)" - WRITE(fileID, cellFormat, advance="no") getCellType(self%cells(e)%obj) + types(e) = getCellType(self%cells(e)%obj) END DO + WRITE(fileID, "(6(I10))") types WRITE(fileID, "(10X, A)") '' WRITE(fileID, "(8X, A)") '' @@ -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)") '' !Write density WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileID,nodeFormat) (output(n)%density, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (output(n)%density, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !Write velocity WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" - WRITE(fileID,nodeFormat) (output(n)%velocity, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (output(n)%velocity, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !Write pressure WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileID,nodeFormat) (output(n)%pressure, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (output(n)%pressure, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !Write temperature WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileID,nodeFormat) (output(n)%temperature, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (output(n)%temperature, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !End node data WRITE(fileID,"(8X,A)") '' @@ -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)") '' 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)") '' - WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(I10))" - WRITE(fileID, cellFormat) (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells) + WRITE(fileID, "(6(I10))") (self%cells(n)%obj%tallyColl(k)%tally(c), n = 1, self%numCells) WRITE(fileID, "(10X, A)") '' 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)") '' !Electric potential WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%phi*Volt_ref, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%phi*Volt_ref, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' !Magnetic Field WRITE(fileID,"(10X,A)") '' - WRITE(nodeFormat, "(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" - WRITE(fileID,nodeFormat) (self%nodes(n)%obj%emData%B*B_ref, n = 1, self%numNodes) + WRITE(fileID,"(6(ES20.6E3))") (self%nodes(n)%obj%emData%B*B_ref, n = 1, self%numNodes) WRITE(fileID,"(10X,A)") '' WRITE(fileID,"(8X,A)") '' !Cell Data + Xi = (/ 0.D0, 0.D0, 0.D0 /) WRITE(fileID,"(8X,A)") '' !Electric field WRITE(fileID,"(10X,A, A, A)") '' - WRITE(cellFormat, "(A,I10, A)") "(", self%numCells, "(3(ES20.6E3)))" - WRITE(fileID, cellFormat) (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) + WRITE(fileID, "(6(ES20.6E3))") (self%cells(n)%obj%gatherElectricField(Xi)*EF_ref, n = 1, self%numCells) WRITE(fileID,"(10X,A)") '' WRITE(fileID,"(8X,A)") '' @@ -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)") '' WRITE(fileIDDeviation,"(10X,A)") '' - WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileIDMean, nodeFormat) (outputMean(n)%density, n = 1, self%numNodes) - WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%density, n = 1, self%numNodes) + 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)") '' WRITE(fileIDDeviation,"(10X,A)") '' !Write velocity WRITE(fileIDMean, "(10X,A)") '' WRITE(fileIDDeviation,"(10X,A)") '' - WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(3(ES20.6E3)))" - WRITE(fileIDMean, nodeFormat) (outputMean(n)%velocity, n = 1, self%numNodes) - WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%velocity, n = 1, self%numNodes) + 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)") '' WRITE(fileIDDeviation,"(10X,A)") '' !Write pressure WRITE(fileIDMean, "(10X,A)") '' WRITE(fileIDDeviation,"(10X,A)") '' - WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileIDMean, nodeFormat) (outputMean(n)%pressure, n = 1, self%numNodes) - WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%pressure, n = 1, self%numNodes) + 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)") '' WRITE(fileIDDeviation,"(10X,A)") '' !Write temperature WRITE(fileIDMean, "(10X,A)") '' WRITE(fileIDDeviation,"(10X,A)") '' - WRITE(nodeFormat,"(A,I10, A)") "(", self%numNodes, "(1(ES20.6E3)))" - WRITE(fileIDMean, nodeFormat) (outputMean(n)%temperature, n = 1, self%numNodes) - WRITE(fileIDDeviation,nodeFormat) (outputDeviation(n)%temperature, n = 1, self%numNodes) + 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)") '' WRITE(fileIDDeviation,"(10X,A)") '' !End node data