diff --git a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90 b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
index f9588a7..3001eac 100644
--- a/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
+++ b/src/modules/mesh/inout/vtu/moduleMeshInputVTU.f90
@@ -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)
 
-        CASE(1)
-          !Edges are nodes, type 1 in VTK
-          self%numEdges = COUNT(types==1)
-        
-        END SELECT
+        self%numEdges = numEdges
 
-        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