Creating of nodes and edges in .vtu format

Moving forward making vtu an independent format.

Now fpakc can generate nodes and edges from vtu input.
Next step is cells.

Some minor corrections in gmsh2 format to unify statements.
The reading of meshes needs a good overhaul.

Testing all geometries with vtu is gonna be fun...

src/modules/init/moduleInput.f90

@@ -860,8 +860,13 @@ MODULE moduleInput
     SUBROUTINE readGeometry(config)
       USE moduleMesh
       USE moduleMeshInputGmsh2, ONLY: initGmsh2
-      USE moduleMeshInputVTU,   ONLY: initVTU, readVTU !TEMPORARY TO TEST VTU OUTPUT
+      USE moduleMeshInputVTU,   ONLY: initVTU
       USE moduleMeshInput0D,    ONLY: init0D
+      USE moduleMesh3DCart
+      USE moduleMesh2DCyl
+      USE moduleMesh2DCart
+      USE moduleMesh1DRad
+      USE moduleMesh1DCart
       USE moduleErrors
       USE moduleOutput
       USE moduleRefParam
@@ -959,7 +964,34 @@ MODULE moduleInput
 
         END SELECT
 
-        !Get the type of mesh
+        !Link the procedure to connect meshes
+        SELECT CASE(mesh%dimen)
+        CASE(3)
+          mesh%connectMesh => connectMesh3DCart
+
+        CASE(2)
+          SELECT CASE(mesh%geometry)
+          CASE("Cyl")
+            mesh%connectMesh => connectMesh2DCyl
+
+          CASE("Cart")
+            mesh%connectMesh => connectMesh2DCart
+
+          END SELECT
+
+        CASE(1)
+          SELECT CASE(mesh%geometry)
+          CASE("Rad")
+            mesh%connectMesh => connectMesh1DRad
+
+          CASE("Cart")
+            mesh%connectMesh => connectMesh1DCart
+
+          END SELECT
+
+        END SELECT
+
+        !Get the format of mesh
         CALL config%get(object // '.meshType', meshFormat, found)
         SELECT CASE(meshFormat)
         CASE ("gmsh2")
@@ -969,14 +1001,22 @@ MODULE moduleInput
 
           END IF
 
+        CASE ("vtu")
+          CALL initVTU(mesh)
+          IF (doubleMesh) THEN
+            CALL initVTU(meshColl)
+
+          END IF
+
+        CASE DEFAULT
+          CALL criticalError('Mesh format ' // meshFormat // ' not defined.', 'readGeometry')
+
         END SELECT
 
         !Reads the mesh file
         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/gmsh2/moduleMeshInputGmsh2.f90

@@ -71,30 +71,26 @@ MODULE moduleMeshInputGmsh2
         SELECT CASE(self%dimen)
         CASE(3)
           ALLOCATE(meshNode3Dcart::self%nodes(n)%obj)
-          self%connectMesh => connectMesh3DCart
 
         CASE(2)
           SELECT CASE(self%geometry)
           CASE("Cyl")
             ALLOCATE(meshNode2DCyl:: self%nodes(n)%obj)
-            self%connectMesh => connectMesh2DCyl
 
           CASE("Cart")
             ALLOCATE(meshNode2DCart:: self%nodes(n)%obj)
-            self%connectMesh => connectMesh2DCart
 
           END SELECT
+
           r(3) = 0.D0
 
         CASE(1)
           SELECT CASE(self%geometry)
           CASE("Rad")
             ALLOCATE(meshNode1DRad:: self%nodes(n)%obj)
-            self%connectMesh => connectMesh1DRad
 
           CASE("Cart")
             ALLOCATE(meshNode1DCart:: self%nodes(n)%obj)
-            self%connectMesh => connectMesh1DCart
 
           END SELECT
           r(2:3) = 0.D0
@@ -111,7 +107,7 @@ MODULE moduleMeshInputGmsh2
       !Reads total number of elements (no nodes)
       READ(10, *) totalNumElem
 
-      !conts edges and volume elements
+      !Count edges and volume elements
       SELECT TYPE(self)
       TYPE IS(meshParticles)
         self%numEdges = 0
@@ -151,7 +147,7 @@ MODULE moduleMeshInputGmsh2
 
       END SELECT
 
-      !Allocates arrays
+      !Allocates array of cells
       ALLOCATE(self%cells(1:self%numCells))
 
       SELECT TYPE(self)
@@ -179,45 +175,30 @@ MODULE moduleMeshInputGmsh2
             END SELECT
 
           CASE (2)
+            ALLOCATE(p(1:2))
+            READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
+            !Associate boundary condition procedure.
+            bt = getBoundaryId(boundaryType)
+
             SELECT CASE(self%geometry)
             CASE("Cyl")
-              ALLOCATE(p(1:2))
-
-              READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
-              !Associate boundary condition procedure.
-              bt = getBoundaryId(boundaryType)
-
               ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
 
             CASE("Cart")
-              ALLOCATE(p(1:2))
-
-              READ(10,*) n, elemType, eTemp, boundaryType, eTemp, p(1:2)
-              !Associate boundary condition procedure.
-              bt = getBoundaryId(boundaryType)
-
               ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)
 
             END SELECT
 
           CASE(1)
+            ALLOCATE(p(1:1))
+            READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
+            !Associate boundary condition 
+            bt = getBoundaryId(boundaryType)
             SELECT CASE(self%geometry)
             CASE("Rad")
-              ALLOCATE(p(1:1))
-
-              READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
-              !Associate boundary condition 
-              bt = getBoundaryId(boundaryType)
-
               ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)
 
             CASE("Cart")
-              ALLOCATE(p(1:1))
-
-              READ(10, *) n, elemType, eTemp, boundaryType, eTemp, p(1)
-              !Associate boundary condition 
-              bt = getBoundaryId(boundaryType)
-
               ALLOCATE(meshEdge1DCart:: self%edges(e)%obj)
 
             END SELECT
@@ -231,9 +212,9 @@ MODULE moduleMeshInputGmsh2
 
       END SELECT
 
-      !Read and initialize volumes
+      !Read and initialize cells
       DO e = 1, self%numCells
-        !Reads the volume according to the geometry
+        !Read the cell according to the geometry
         SELECT CASE(self%dimen)
         CASE(3)
           READ(10, *) n, elemType

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

@@ -124,7 +124,6 @@ MODULE moduleMeshInputVTU
       DO WHILE (iStart < nData)
         iStart = iStart + 1
         iEnd   = iStart - 1 + block
-        PRINT *, iStart, iEnd
         IF (iEnd > nData) THEN
           iEnd = nData
 
@@ -173,6 +172,9 @@ MODULE moduleMeshInputVTU
       INTEGER:: numNodes, numElements
       INTEGER, ALLOCATABLE, DIMENSION(:):: entitiesID, offsets, connectivity, types
       REAL(8), ALLOCATABLE, DIMENSION(:):: coordinates
+      INTEGER:: n, e, c
+      INTEGER, ALLOCATABLE:: p(:)
+      INTEGER:: bt
 
       fileID = 10
 
@@ -200,6 +202,8 @@ MODULE moduleMeshInputVTU
       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
@@ -217,7 +221,164 @@ MODULE moduleMeshInputVTU
       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)
+
+        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
+        CALL self%nodes(n)%obj%init(n, r)
+
+      END DO
+
+      !Count the number of edges
+      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)
+
+        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
+      e = 0
+      c = 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(meshEdge3DCartTria:: self%edges(e)%obj)
+              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
+              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))
+
+              SELECT CASE(self%geometry)
+              CASE("Cyl")
+                ALLOCATE(meshEdge2DCyl:: self%edges(e)%obj)
+
+              CASE("Cart")
+                ALLOCATE(meshEdge2DCart:: self%edges(e)%obj)
+
+              END SELECT
+
+              !Allocate edge
+              CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
+              DEALLOCATE(p)
+
+            END IF
+
+          CASE(1)
+            IF (types(n) == 3) THEN
+              e = e + 1
+              ALLOCATE(p(1:1))
+              p(1) = connectivity(offsets(n))
+
+              !Associate boundary condition procedure.
+              bt = getBoundaryId(entitiesID(n))
+              SELECT CASE(self%geometry)
+              CASE("Rad")
+                ALLOCATE(meshEdge1DRad:: self%edges(e)%obj)
+
+              CASE("Cart")
+                ALLOCATE(meshEdge1DCart:: self%edges(e)%obj)
+
+              END SELECT
+
+              CALL self%edges(e)%obj%init(n, p, bt, entitiesID(n))
+              DEALLOCATE(p)
+
+            END IF
+
+          END SELECT
+
+        END DO
+
+      END SELECT
+        
     END SUBROUTINE readVTU
 
     SUBROUTINE readInitialVTU(filename, density, velocity, temperature)

src/modules/solver/moduleSolver.f90

@@ -514,7 +514,6 @@ MODULE moduleSolver
       USE moduleSpecies
       USE moduleCompTime
       USE moduleProbe
-      USE moduleMeshOutputVTU !TEMPORARY TO TEST VTU OUTPUT
       IMPLICIT NONE
 
       INTEGER, INTENT(in):: t
@@ -528,11 +527,8 @@ MODULE moduleSolver
 
         CALL outputProbes(t)
         CALL mesh%printOutput(t)
-        CALL printOutputVTU(mesh, t) !TEMPORARY TO TEST VTU OUTPUT
         IF (ASSOCIATED(meshForMCC)) CALL meshForMCC%printColl(t)
-        IF (ASSOCIATED(meshForMCC)) CALL printCollVTU(meshForMCC,t) !TEMPORARY TO TEST VTU OUTPUT
         CALL mesh%printEM(t)
-        CALL printEMVTU(mesh, t) !TEMPORARY TO TEST VTU OUTPUT
         WRITE(*, "(5X,A21,I10,A1,I10)") "t/tFinal: ", t, "/", tFinal
         WRITE(*, "(5X,A21,I10)")        "Particles: ", nPartOld
         IF (t == 0) THEN
@@ -565,7 +561,6 @@ MODULE moduleSolver
       !Output average values
       IF (useAverage .AND. t == tFinal) THEN
         CALL mesh%printAverage()
-        CALL printAverageVTU(mesh) !TEMPORARY TO TEST VTU OUTPUT
 
       END IF