Improvement to vtu read subroutine

Some improvements to reduce code repetition when reading vtu mesh.
This commit is contained in:
Jorge Gonzalez 2023-02-06 19:18:04 +01:00
commit a22099ee87

View file

@ -169,7 +169,7 @@ MODULE moduleMeshInputVTU
REAL(8):: r(1:3) !3 generic coordinates
INTEGER:: fileID, error, found
CHARACTER(LEN=256):: line
INTEGER:: numNodes, numElements
INTEGER:: numNodes, numElements, numEdges
INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types
REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates
INTEGER:: n, e, c
@ -223,6 +223,7 @@ MODULE moduleMeshInputVTU
!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))
@ -238,6 +239,7 @@ MODULE moduleMeshInputVTU
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)
@ -266,41 +268,44 @@ MODULE moduleMeshInputVTU
r(2:3) = 0.D0
END SELECT
!Init node
CALL self%nodes(n)%obj%init(n, r)
END DO
!Count the number of edges
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)
SELECT CASE(self%dimen)
CASE(3)
!Edges are triangles, type 5 in VTK
self%numEdges = COUNT(types==5)
CASE(2)
!Edges are segments, type 3 in VTK
self%numEdges = COUNT(types==3)
self%numEdges = numEdges
CASE(1)
!Edges are nodes, type 1 in VTK
self%numEdges = COUNT(types==1)
END SELECT
self%numCells = numElements - self%numEdges
!Allocate array of edges
ALLOCATE(self%edges(1:self%numEdges))
TYPE IS(meshCollisions)
self%numCells = numElements
END SELECT
!Allocates array of cells
ALLOCATE(self%cells(1:self%numCells))
!Read edges
!Read edges only if mesh is for tracking particles, not for collisions
e = 0
SELECT TYPE(self)
TYPE IS(meshParticles)
@ -309,9 +314,8 @@ MODULE moduleMeshInputVTU
CASE(3)
IF (types(n) == 5) THEN
e = e + 1
ALLOCATE(meshEdge3DCartTria:: self%edges(e)%obj)
ALLOCATE(p(1:3))
ALLOCATE(p(1:3))
p(1) = connectivity(offsets(n) - 2)
p(2) = connectivity(offsets(n) - 1)
p(3) = connectivity(offsets(n))
@ -320,6 +324,9 @@ MODULE moduleMeshInputVTU
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)
@ -335,6 +342,7 @@ MODULE moduleMeshInputVTU
!Associate boundary condition procedure.
bt = getBoundaryId(entitiesID(n))
!Allocate edge
SELECT CASE(self%geometry)
CASE("Cyl")
ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
@ -344,7 +352,7 @@ MODULE moduleMeshInputVTU
END SELECT
!Allocate edge
!Init edge
CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
DEALLOCATE(p)
@ -358,6 +366,8 @@ MODULE moduleMeshInputVTU
!Associate boundary condition procedure.
bt = getBoundaryId(entitiesID(n))
!Allocate edge
SELECT CASE(self%geometry)
CASE("Rad")
ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)
@ -367,6 +377,7 @@ MODULE moduleMeshInputVTU
END SELECT
!Init edge
CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
DEALLOCATE(p)
@ -387,112 +398,98 @@ MODULE moduleMeshInputVTU
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(self%geometry)
CASE("Cyl")
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))
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)
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(meshCell2DCylQuad:: self%cells(c)%obj)
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
CASE("Cart")
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))
CASE("Cart")
ALLOCATE(meshCell2DCartTria:: self%cells(c)%obj)
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
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))
!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)
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
!Init cell
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
CASE(1)
SELECT CASE(self%geometry)
CASE("Rad")
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))
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)
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
CASE("Cart")
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))
CASE("Cart")
ALLOCATE(meshCell1DCartSegm:: self%cells(c)%obj)
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
!Init cell
CALL self%cells(c)%obj%init(c, p, self%nodes)
DEALLOCATE(p)
END SELECT
END SELECT