diff --git a/anisotropy/meshing/mesh.py b/anisotropy/meshing/mesh.py
index 11cc720..c374528 100644
--- a/anisotropy/meshing/mesh.py
+++ b/anisotropy/meshing/mesh.py
@@ -6,6 +6,8 @@ from netgen.occ import OCCGeometry
 from netgen import meshing
 from numpy import array
 import os
+from .utils import extractPoints, extractCells
+import meshio
 
 
 class NoGeometrySpecified(Exception):
@@ -114,21 +116,55 @@ class Mesh(object):
         
         return out, err, returncode
 
-    def doubleExport(self):
-        pass
+    def to_meshio(self): 
+        points = extractPoints(self.mesh.Points())
+        cells = []
+
+        if len(self.mesh.Elements1D()) > 0:
+            cells.extend([ cells_ for cells_ in extractCells(1, self.mesh.Elements1D()).items() ])   
+
+        if len(self.mesh.Elements2D()) > 0:
+            cells.extend([ cells_ for cells_ in extractCells(2, self.mesh.Elements2D()).items() ])   
+
+        if len(self.mesh.Elements3D()) > 0:
+            cells.extend([ cells_ for cells_ in extractCells(3, self.mesh.Elements3D()).items() ])   
+
+        return meshio.Mesh(points, cells)
+
+    @staticmethod
+    def volumeTetra(points: array) -> float:
+        return 1 / 6 * linalg.det(numpy.append(points.transpose(), numpy.array([[1, 1, 1, 1]]), axis = 0))
 
     @property
     def volumes(self) -> array:
-        return array(self.mesh.Elements3D())
+        points = []
+
+        for cells in self.cells:
+            if cells.dim == 3:
+                points.extend([ self.mesh.points[cell] for cell in cells.data ])
+
+        return array(points)
 
     @property
     def faces(self) -> array:
-        return array(self.mesh.Elements2D())
+        points = []
+
+        for cells in self.cells:
+            if cells.dim == 2:
+                points.extend([ self.mesh.points[cell] for cell in cells.data ])
+
+        return array(points)
 
     @property
     def edges(self) -> array:
-        # NOTE: returns zero elements for imported mesh 
-        return array(self.mesh.Elements1D())
+        points = []
+
+        for cells in self.cells:
+            if cells.dim == 1:
+                points.extend([ self.mesh.points[cell] for cell in cells.data ])
+
+        return array(points)
+
 
 
 # tetras = numpy.array([ [ [ vertex for vertex in mesh[index] ] for index in element.vertices ] for element in self.mesh.Elements3D() ])
diff --git a/anisotropy/meshing/utils.py b/anisotropy/meshing/utils.py
new file mode 100644
index 0000000..8f0e619
--- /dev/null
+++ b/anisotropy/meshing/utils.py
@@ -0,0 +1,31 @@
+from meshio._mesh import topological_dimension
+from meshio._common import num_nodes_per_cell
+from numpy import array
+
+
+def detectTopology(dimension: dict, num_nodes: dict):
+    for dim in topological_dimension.keys():
+        for num in num_nodes_per_cell.keys():
+            if topological_dimension[dim] == dimension and num_nodes_per_cell[num] == num_nodes and dim == num:
+                return dim
+
+
+def extractPoints(points):
+    return array([ point.p for point in points ], dtype = float)
+
+
+def extractCells(dimension: int, elements):
+    cellsNew = {}
+
+    for cell in elements:
+        cellTopo = detectTopology(dimension, len(cell.points))
+        #   shift indicies, they should starts from zero
+        cellNew = array([ pointId.nr for pointId in cell.points ], dtype = int) - 1
+
+        if cellsNew.get(cellTopo):
+            cellsNew[cellTopo].append(cellNew)
+
+        else:
+            cellsNew[cellTopo] = [ cellNew ] 
+        
+    return cellsNew
diff --git a/playground/meshConvert.py b/playground/meshConvert.py
index d595cf8..67b4ff5 100644
--- a/playground/meshConvert.py
+++ b/playground/meshConvert.py
@@ -1,33 +1,47 @@
 from netgen.meshing import Mesh as netgenMesh
 from meshio import Mesh as meshioMesh
+from meshio._mesh import topological_dimension
+from meshio._common import num_nodes_per_cell
+from numpy import array
 
 meshfile = "mesh.mesh"
 
 mesh = netgenMesh()
 mesh.Load(meshfile)
 
-topology3d = {
-    4: "tetra"
-}
+def detectTopology(dimension: dict, num_nodes: dict):
+    for dim in topological_dimension.keys():
+        for num in num_nodes_per_cell.keys():
+            if topological_dimension[dim] == dimension and num_nodes_per_cell[num] == num_nodes and dim == num:
+                return dim
 
-pointsNew = []
-cellsNew = {}
+def extractCells(dimension: int, elements):
+    cellsNew = {}
 
-for point in mesh.Points():
-    pointsNew.append(list(point.p))
+    for cell in elements:
+        cellTopo = detectTopology(dimension, len(cell.points))
+        #   shift indicies, they should starts from zero
+        cellNew = array([ pointId.nr for pointId in cell.points ], dtype = int) - 1
 
-for cell in mesh.Elements3D():
-    cellTopo = topology3d[len(cell.points)]
-    cellNew = []
+        if cellsNew.get(cellTopo):
+            cellsNew[cellTopo].append(cellNew)
 
-    for pointId in cell.points:
-        cellNew.append(pointId.nr)
+        else:
+            cellsNew[cellTopo] = [ cellNew ] 
+        
+    return cellsNew
 
-    if cellsNew.get(cellTopo):
-        cellsNew[cellTopo].append(cellNew)
+def extractPoints(points):
+    return array([ point.p for point in mesh.Points() ], dtype = float)
 
-    else:
-        cellsNew[cellTopo] = [ cellNew ]    
+points = extractPoints(mesh.Points())
+cells1d = extractCells(1, mesh.Elements1D())   
+cells2d = extractCells(2, mesh.Elements2D())   
+cells3d = extractCells(3, mesh.Elements3D())   
+cells = [
+    *[ e for e in cells1d.items() ],
+    *[ e for e in cells2d.items() ],
+    *[ e for e in cells3d.items() ]
+]
 
-cellsMeshio = [ cell for cell in cellsNew.items() ]
-meshNew = meshioMesh(pointsNew, cellsMeshio)
\ No newline at end of file
+meshNew = meshioMesh(points, cells)
\ No newline at end of file