diff --git a/examples/dem2.py b/examples/dem2.py
new file mode 100644
index 0000000000000000000000000000000000000000..7ebb47b1a335c4e0728a62500ecb0e2adcd81402
--- /dev/null
+++ b/examples/dem2.py
@@ -0,0 +1,111 @@
+import pairs
+import sys
+
+
+def linear_spring_dashpot(i, j):
+    penetration_depth = squared_distance(i, j) - radius[i] - radius[j]
+    skip_when(penetration_depth >= 0.0)
+
+    delta_ij = delta(i, j)
+    contact_normal = delta_ij / length(delta_ij)
+    k = radius[j] + 0.5 * penetration_depth
+    contact_point = position[j] + contact_normal * k
+
+    rel_vel = -velocity_wf[i] - velocity_wf[j]
+    rel_vel_n = dot(rel_vel, contact_normal) * contact_normal
+    rel_vel_t = rel_vel - rel_vel_n
+
+    fN = stiffness_norm[i, j] * (-penetration_depth) * contact_normal + damping_norm[i, j] * rel_vel_n;
+
+    tan_spring_disp = tangential_spring_displacement[i, j]
+    impact_vel_magnitude = impact_velocity_magnitude[i, j]
+    impact_magnitude = select(impact_vel_magnitude > 0.0, impact_vel_magnitude, length(rel_vel))
+    sticking = is_sticking[i, j]
+
+    rotated_tan_disp = tan_spring_disp - contact_normal * (contact_normal * tan_spring_disp)
+    new_tan_spring_disp = dt * rel_vel_t + \
+                          select(squared_length(rotated_tan_disp) <= 0.0,
+                                 zero_vector(),
+                                 rotated_tan_disp * length(tan_spring_disp) / length(rotated_tan_disp))
+
+    fTLS = stiffness_tan[i, j] * new_tan_spring_disp + damping_tan[i, j] * rel_vel_t
+    fTLS_len = length(fTLS)
+    fTLS_inv_len = 1.0 / fTLS_len
+    t = select(fTLS_len > 0, fTLS / fTLS_inv_len, zero_vector())
+
+    f_friction_abs_static = friction_static[i, j] * length(fN)
+    f_friction_abs_dynamic = friction_dynamic[i, j] * length(fN)
+    tan_vel_threshold = 1e-8
+
+    cond1 = sticking == 1 and length(rel_vel_t) < tan_vel_threshold and fTLS_len < f_friction_abs_static
+    cond2 = sticking == 1 and fTLS_len < f_friction_abs_dynamic
+    f_friction_abs = select(cond1, f_friction_abs_static, f_friction_abs_dynamic)
+    n_sticking = select(cond1 or cond2 or fTLS_len < f_friction_abs_dynamic, 1, 0)
+
+    if not cond1 and not cond2 and stiffness_tan[i, j] > 0.0:
+        tangential_spring_displacement[i, j] = \
+            (f_friction_abs * t - damping_tan[i, j] * rel_vel_t) / stiffness_tan[i, j]
+
+    else:
+        tangential_spring_displacement[i, j] = new_tan_spring_disp
+
+    impact_velocity_magnitude[i, j] = impact_magnitude
+    is_sticking[i, j] = n_sticking
+
+    fTabs = min(fTLS_len, f_friction_abs)
+    fT = fTabs * t
+    force[i] += fN + fT
+
+
+def euler(i):
+    velocity[i] += dt * force[i] / mass[i]
+    position[i] += dt * velocity[i]
+
+
+cmd = sys.argv[0]
+target = sys.argv[1] if len(sys.argv[1]) > 1 else "none"
+if target != 'cpu' and target != 'gpu':
+    print(f"Invalid target, use {cmd} <cpu/gpu>")
+
+
+dt = 0.005
+cutoff_radius = 2.5
+skin = 0.3
+ntypes = 4
+stiffness_norm = 1.0
+stiffness_tan = 1.0
+damping_norm = 1.0
+damping_tan = 1.0
+friction_static = 1.0
+friction_dynamic = 1.0
+
+psim = pairs.simulation("dem", debug=True)
+psim.add_position('position')
+psim.add_property('mass', pairs.double(), 1.0)
+psim.add_property('velocity', pairs.vector())
+psim.add_property('velocity_wf', pairs.vector())
+psim.add_property('force', pairs.vector(), vol=True)
+psim.add_property('radius', pairs.double(), 1.0)
+psim.add_feature('type', ntypes)
+psim.add_feature_property('type', 'stiffness_norm', pairs.double(), [stiffness_norm for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'stiffness_tan', pairs.double(), [stiffness_tan for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'damping_norm', pairs.double(), [damping_norm for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'damping_tan', pairs.double(), [damping_tan for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'friction_static', pairs.double(), [friction_static for i in range(ntypes * ntypes)])
+psim.add_feature_property('type', 'friction_dynamic', pairs.double(), [friction_dynamic for i in range(ntypes * ntypes)])
+psim.add_contact_property('is_sticking', pairs.int32(), False)
+psim.add_contact_property('tangential_spring_displacement', pairs.vector(), [0.0, 0.0, 0.0])
+psim.add_contact_property('impact_velocity_magnitude', pairs.double(), 0.0)
+
+psim.read_particle_data("data/fluidized_bed.input", ['mass', 'position', 'velocity'])
+psim.build_neighbor_lists(cutoff_radius + skin)
+psim.vtk_output(f"output/test_{target}")
+psim.compute(linear_spring_dashpot, cutoff_radius, symbols={'dt': dt})
+psim.compute(euler, symbols={'dt': dt})
+
+if target == 'gpu':
+    psim.target(pairs.target_gpu())
+else:
+    psim.target(pairs.target_cpu())
+
+psim.generate()
diff --git a/src/pairs/graph/graphviz.py b/src/pairs/graph/graphviz.py
index ba2836003c8bd9de5e362b26c70423b0b3ceee08..08210da886c9b5b3ca7bdc49d661a6ad891da757 100644
--- a/src/pairs/graph/graphviz.py
+++ b/src/pairs/graph/graphviz.py
@@ -1,43 +1,54 @@
 from graphviz import Digraph
 from pairs.ir.arrays import Array
 from pairs.ir.bin_op import BinOp, Decl
+from pairs.ir.features import Feature, FeatureProperty
 from pairs.ir.lit import Lit
 from pairs.ir.loops import Iter
-from pairs.ir.properties import Property
+from pairs.ir.properties import Property, ContactProperty
 from pairs.ir.variables import Var
 from pairs.ir.visitor import Visitor
 
 
 class ASTGraph:
+    last_reference = 0
+
     def __init__(self, ast_node, filename, ref="AST", max_depth=0):
         self.graph = Digraph(ref, filename=filename, node_attr={'color': 'lightblue2', 'style': 'filled'})
         self.graph.attr(size='6,6')
-        self.visitor = Visitor(ast_node, max_depth=max_depth)
+        self.ast = ast_node
 
-    def render(self):
-        def generate_edges_for_node(ast_node, graph, generated):
-            node_id = id(ast_node)
-            if not isinstance(ast_node, Decl) and node_id not in generated:
-                node_ref = f"n{id(ast_node)}"
-                generated.append(node_id)
-                graph.node(node_ref, label=ASTGraph.get_node_label(ast_node))
-                for child in ast_node.children():
-                    if not isinstance(child, Decl):
-                        child_ref = f"n{id(child)}"
-                        graph.node(child_ref, label=ASTGraph.get_node_label(child))
-                        graph.edge(node_ref, child_ref)
-
-        generated = []
-        for node in self.visitor:
-            generate_edges_for_node(node, self.graph, generated)
+    def new_unique_reference():
+        ASTGraph.last_reference += 1
+        return f"unique_ref{ASTGraph.last_reference}"
+
+    def generate_edges_for_node(ast_node, graph, generated):
+        node_ref = ASTGraph.get_node_reference(ast_node)
+        if not isinstance(ast_node, Decl) and node_ref not in generated:
+            generated.append(node_ref)
+            graph.node(node_ref, label=ASTGraph.get_node_label(ast_node))
+
+            for child in ast_node.children():
+                if not isinstance(child, Decl):
+                    child_ref = ASTGraph.generate_edges_for_node(child, graph, generated)
+                    graph.edge(node_ref, child_ref)
 
+        return node_ref
+
+    def render(self):
+        ASTGraph.generate_edges_for_node(self.ast, self.graph, [])
         self.graph.render()
 
     def view(self):
         self.graph.view()
 
+    def get_node_reference(ast_node):
+        if isinstance(ast_node, (Array, Property, Var, ContactProperty, FeatureProperty, Iter)):
+            return ASTGraph.new_unique_reference()
+
+        return f"n{id(ast_node)}"
+
     def get_node_label(ast_node):
-        if isinstance(ast_node, (Array, Property, Var)):
+        if isinstance(ast_node, (Array, Property, Var, ContactProperty, FeatureProperty)):
             return ast_node.name()
 
         if isinstance(ast_node, BinOp):
diff --git a/src/pairs/ir/mutator.py b/src/pairs/ir/mutator.py
index ff59c47911ed55a1374398043d45352033e126e3..416e77b4caff0f1d45f0df50a51492bb067503b6 100644
--- a/src/pairs/ir/mutator.py
+++ b/src/pairs/ir/mutator.py
@@ -2,6 +2,7 @@ class Mutator:
     def __init__(self, ast=None, max_depth=0):
         self.ast = ast
         self.max_depth = 0
+        self.mutated_vector_expressions = []
 
     def set_ast(self, ast):
         self.ast = ast
@@ -180,7 +181,13 @@ class Mutator:
         return ast_node
 
     def mutate_VectorAccess(self, ast_node):
-        ast_node.expr = self.mutate(ast_node.expr)
+        # Traversing the expressions for all vector accesses makes the compilation
+        # significantly slower, when it is necessary, a specific method can be used
+        # for the transformation
+        if id(ast_node.expr) not in self.mutated_vector_expressions:
+            ast_node.expr = self.mutate(ast_node.expr)
+            self.mutated_vector_expressions.append(id(ast_node.expr))
+
         return ast_node
 
     def mutate_While(self, ast_node):
diff --git a/src/pairs/ir/visitor.py b/src/pairs/ir/visitor.py
index 2aa56c874e9d34875aab24072562c526c246bed1..14124aed2e7c8452772ee9c9b4c027ae3392282c 100644
--- a/src/pairs/ir/visitor.py
+++ b/src/pairs/ir/visitor.py
@@ -6,6 +6,7 @@ class Visitor:
         self.ast = ast
         self.max_depth = max_depth
         self.breadth_first = breadth_first
+        self.visited_vector_expressions = []
 
     def set_ast(self, ast):
         self.ast = ast
@@ -35,6 +36,15 @@ class Visitor:
                 if method is None:
                     self.visit(node.children())
 
+    #def visit_VectorAccess(self, ast_node):
+        # Traversing the expressions for all vector accesses makes the compilation
+        # significantly slower, when it is necessary, a specific method can be used
+        # for the analysis
+    #    expr_id = id(ast_node.expr)
+    #    if expr_id not in self.visited_vector_expressions:
+    #        self.visit(ast_node.children())
+    #        self.visited_vector_expressions.append(expr_id)
+
     def visit_children(self, ast_node):
         self.visit(ast_node.children())
 
diff --git a/src/pairs/mapping/keywords.py b/src/pairs/mapping/keywords.py
index 36e83919c9d979ac7ef2209838667f20021385e5..cbcd3739a5cf748dd82ffdf680d075748e2c3393 100644
--- a/src/pairs/mapping/keywords.py
+++ b/src/pairs/mapping/keywords.py
@@ -64,7 +64,7 @@ class Keywords:
     def keyword_normalized(self, args):
         assert len(args) == 1, "normalized() keyword requires one parameter!"
         vector = args[0]
-        assert vector.type() == Types.Vector, "length(): Argument must be a vector!"
+        assert vector.type() == Types.Vector, "normalized(): Argument must be a vector!"
         length = self.keyword_length([vector])
         inv_length = Lit(self.sim, 1.0) / length
         return Select(self.sim, length > Lit(self.sim, 0.0), vector * inv_length, ZeroVector(self.sim))
diff --git a/src/pairs/sim/simulation.py b/src/pairs/sim/simulation.py
index c28d54c4aac3a953237021046c71558ee51c80b8..427462c9499f8f64e1a2fc76d38342ab27a0e115 100644
--- a/src/pairs/sim/simulation.py
+++ b/src/pairs/sim/simulation.py
@@ -10,7 +10,7 @@ from pairs.ir.properties import Properties, ContactProperties
 from pairs.ir.symbols import Symbol
 from pairs.ir.types import Types
 from pairs.ir.variables import Variables
-#from pairs.graph.graphviz import ASTGraph
+from pairs.graph.graphviz import ASTGraph
 from pairs.mapping.funcs import compute
 from pairs.sim.arrays import ArraysDecl
 from pairs.sim.cell_lists import CellLists, CellListsBuild, CellListsStencilBuild