From e3eeb97f4863117a8d138c0ca9b2843e5f7cfdda Mon Sep 17 00:00:00 2001
From: Jorge Gonzalez <jorge.gonzalez@upm.es>
Date: Sun, 5 Feb 2023 21:33:03 +0100
Subject: [PATCH] First working version!

First complete implementation of .vtu format.
Still a lot of things to improve but right now fpakc can read a vtu mesh
and write the output in vtu.

Still to test:
  Multiple geometries.
  Double mesh.
---
 .../mesh/inout/vtu/moduleMeshInputVTU.f90     | 125 +++++++++++++++++-
 1 file changed, 124 insertions(+), 1 deletion(-)

diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
index 3a72d12..f9588a7 100644
--- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
@@ -302,7 +302,6 @@ MODULE moduleMeshInputVTU
 
       !Read edges
       e = 0
-      c = 0
       SELECT TYPE(self)
       TYPE IS(meshParticles)
         DO n = 1, numElements
@@ -378,7 +377,131 @@ MODULE moduleMeshInputVTU
         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(meshCell3DCartTetra::  self%cells(c)%obj)
+
+            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))
+              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))
+              ALLOCATE(meshCell2DCartTria:: 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(meshCell2DCartQuad:: self%cells(c)%obj)
+
+              CALL self%cells(c)%obj%init(c, p, self%nodes)
+              DEALLOCATE(p)
+
+            END SELECT
+
+          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))
+              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))
+              ALLOCATE(meshCell1DCartSegm:: self%cells(c)%obj)
+
+              CALL self%cells(c)%obj%init(c, p, self%nodes)
+              DEALLOCATE(p)
+
+            END SELECT
+
+          END SELECT
+
+        END SELECT
+
+      END DO
         
+      !Call mesh connectivity
+      CALL self%connectMesh
+
     END SUBROUTINE readVTU
 
     SUBROUTINE readInitialVTU(filename, density, velocity, temperature)