JGonzalez
1b1a574edc
Now we can read an initial condition as we do with other formats. Only documentation left.
555 lines
15 KiB
Fortran
555 lines
15 KiB
Fortran
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
|
|
USE moduleBoundary
|
|
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(:)
|
|
INTEGER:: bt
|
|
|
|
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))
|
|
|
|
!Associate boundary condition procedure.
|
|
bt = getBoundaryId(entitiesID(n))
|
|
|
|
!Allocate edge
|
|
ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj)
|
|
|
|
!Init edge
|
|
CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
|
|
DEALLOCATE(p)
|
|
|
|
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))
|
|
|
|
!Associate boundary condition procedure.
|
|
bt = getBoundaryId(entitiesID(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
|
|
|
|
!Init edge
|
|
CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
|
|
DEALLOCATE(p)
|
|
|
|
END IF
|
|
|
|
CASE(1)
|
|
IF (types(n) == 1) THEN
|
|
e = e + 1
|
|
ALLOCATE(p(1:1))
|
|
p(1) = connectivity(offsets(n))
|
|
|
|
!Associate boundary condition procedure.
|
|
bt = getBoundaryId(entitiesID(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
|
|
|
|
!Init edge
|
|
CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
|
|
DEALLOCATE(p)
|
|
|
|
END IF
|
|
|
|
END SELECT
|
|
|
|
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
|