diff --git a/Project.toml b/Project.toml
index cf0e42c..1f39869 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FieldTracer"
 uuid = "05065131-34d1-456a-ba22-faf972bb2934"
 authors = ["Hongyang Zhou <hyzhou@umich.edu>"]
-version = "0.1.1"
+version = "0.1.2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -10,7 +10,7 @@ PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 
 [compat]
-Meshes = "0.12"
+Meshes = "0.13"
 PyPlot = "2.9"
 Requires = "1.1"
 julia = "1"
diff --git a/src/unstructured2d.jl b/src/unstructured2d.jl
index e6bb4db..5119b41 100644
--- a/src/unstructured2d.jl
+++ b/src/unstructured2d.jl
@@ -24,7 +24,7 @@ function trace(mesh::SimpleMesh, vx, vy, xstart, ystart; maxIter=1000, maxLen=10
 
       for it = 1:maxIter
          Pfar = Pnow + Vec2f(vx[cellID]*Δ, vy[cellID]*Δ)
-         nodes = mesh.connec[cellID].list
+         element = getelement(mesh, cellID)
          
          if mesh[cellID] isa Quadrangle
             nEdge = 4
@@ -36,13 +36,13 @@ function trace(mesh::SimpleMesh, vx, vy, xstart, ystart; maxIter=1000, maxLen=10
          # Loop over edges
          for i = 1:nEdge
             j = i % nEdge + 1 
-            P = Segment(mesh.points[nodes[i]], mesh.points[nodes[j]]) ∩ ray
+            P = Segment(element.vertices[i], element.vertices[j]) ∩ ray
             if P isa Point
                P⁺ = P + Vec2f(vx[cellID]*ϵ, vy[cellID]*ϵ)
                cellIDNew = getCellID(mesh, P⁺)
                break
             elseif P isa Segment
-               P⁺ = mesh.points[nodes[j]] + Vec2f(vx[cellID]*ϵ, vy[cellID]*ϵ)
+               P⁺ = element.vertices[j] + Vec2f(vx[cellID]*ϵ, vy[cellID]*ϵ)
                cellIDNew = getCellID(mesh, P⁺)
                break
             end
@@ -70,25 +70,21 @@ function trace(mesh::SimpleMesh, vx, vy, xstart, ystart; maxIter=1000, maxLen=10
    xStream, yStream
 end
 
-"""
-    getCellID(mesh::SimpleMesh, x, y)
-
-Return cell ID on the unstructured mesh.
-"""
+"Return cell ID on the unstructured mesh."
 function getCellID(mesh::SimpleMesh, point::Point2)
-   for i = 1:length(mesh.connec)
-      nodes = mesh.connec[i].list
-      if mesh[i] isa Triangle
-         if point ∈ Triangle(mesh.points[nodes[1]], mesh.points[nodes[2]],
-            mesh.points[nodes[3]])
-            return i
-         end
-      else # Quadrangle
-         if point ∈ Quadrangle(mesh.points[nodes[1]], mesh.points[nodes[2]],
-            mesh.points[nodes[3]], mesh.points[nodes[4]])
-            return i
-         end
+   for (i, element) in enumerate(elements(mesh))
+      if point ∈ element
+         return i
       end
    end
    return 0 # out of mesh boundary
+end
+
+"Return the `cellID`th element of the mesh."
+function getelement(mesh, cellID)
+   for (i, element) in enumerate(elements(mesh))
+      if i == cellID
+         return element
+      end
+   end
 end
\ No newline at end of file