diff --git a/geosolver/clsolver.py b/geosolver/clsolver.py
index af1b090..99dcc52 100644
--- a/geosolver/clsolver.py
+++ b/geosolver/clsolver.py
@@ -927,5 +927,18 @@ class ClusterSolver(Notifier):
     ##                raise "cluster determined by more than one method"
     ##    return result
 
+
+    ##def _all_sources_constraint_in_cluster(self, constraint, cluster):
+    ##    if not self._contains_constraint(cluster, constraint):
+    ##        return Set()
+    ##    elif self._is_atomic(cluster):
+    ##        return Set([cluster])
+    ##    else:
+    ##        method = self._determining_method(cluster)
+    ##        sources = Set()
+    ##        for inp in method.input_clusters():
+    ##            sources.union_update(self._all_sources_constraint_in_cluster(constraint, inp))
+    ##        return sources
+     
 # class ClusterSolver
 
diff --git a/geosolver/clsolver3D.py b/geosolver/clsolver3D.py
index b38a438..b085ead 100644
--- a/geosolver/clsolver3D.py
+++ b/geosolver/clsolver3D.py
@@ -88,6 +88,7 @@ class ClusterSolver3D(ClusterSolver):
 
     # overriding ClusterSolver.set_root
     def set_root(self, cluster):
+        """Set root cluster, used for positionig and orienting the solutions"""
         diag_print("set root "+str(self.rootcluster), "clsolver3D")
         if self.rootcluster != None:
             oldrootvar = rootname(self.rootcluster)
@@ -99,18 +100,7 @@ class ClusterSolver3D(ClusterSolver):
 
     # ------------ INTERNALLY USED METHODS --------
 
-    def _all_sources_constraint_in_cluster(self, constraint, cluster):
-        if not self._contains_constraint(cluster, constraint):
-            return Set()
-        elif self._is_atomic(cluster):
-            return Set([cluster])
-        else:
-            method = self._determining_method(cluster)
-            sources = Set()
-            for inp in method.input_clusters():
-                sources.union_update(self._all_sources_constraint_in_cluster(constraint, inp))
-            return sources
-     
+    
     # --------------
     # search methods
     # --------------
@@ -212,15 +202,19 @@ class ClusterSolver3D(ClusterSolver):
         self._add_cluster(output)
         self._add_method(merge)
         # remove input clusters from top_level
-        if not (hasattr(merge,"noremove") and merge.noremove == True):
-            merge.restore_toplevel = []    # make restore list in method
-            for cluster in merge.input_clusters():
-                if num_constraints(cluster.intersection(output)) >= num_constraints(cluster): 
-                    diag_print("remove from top-level: "+str(cluster),"clsolver3D")
-                    self._rem_top_level(cluster) 
-                    merge.restore_toplevel.append(cluster)
-                else:
-                    diag_print("keep top-level: "+str(cluster),"clsolver3D")
+        merge.restore_toplevel = []    # make restore list in method
+        for cluster in merge.input_clusters():
+            # do not remove rigids from toplevel if method does not consider root
+            if isinstance(cluster, Rigid):
+                if hasattr(merge,"noremove") and merge.noremove == True:
+                    continue
+            # remove input clusters when all its constraints are in output cluster 
+            if num_constraints(cluster.intersection(output)) >= num_constraints(cluster): 
+                diag_print("remove from top-level: "+str(cluster),"clsolver3D")
+                self._rem_top_level(cluster) 
+                merge.restore_toplevel.append(cluster)
+            else:
+                diag_print("keep top-level: "+str(cluster),"clsolver3D")
         # add method to determine root-variable
         self._add_root_method(merge.input_clusters(),merge.outputs()[0])
         # add solution selection methods
@@ -544,7 +538,7 @@ class DeriveDAD(ClusterMethod):
         return solutions
 
 class DeriveADD(ClusterMethod):
-    """Represents a merging of one distance and to distances"""
+    """Represents a merging of one angle and two distances"""
     def __init__(self, map):
         # check inputs
         self.a_cab = map["$a_cab"]
@@ -636,7 +630,7 @@ class DeriveAA(ClusterMethod):
         return solutions
 
 class MergeSR(ClusterMethod):
-    """Merge a Rigid from a Scalabe and a Rigid sharing two points"""
+    """Merge a Scalabe and a Rigid sharing two points"""
     def __init__(self, map):
         # check inputs
         in1 = map["$r"]
diff --git a/geosolver/cluster.py b/geosolver/cluster.py
index 9c1fb5f..156ad29 100644
--- a/geosolver/cluster.py
+++ b/geosolver/cluster.py
@@ -7,6 +7,7 @@ from multimethod import MultiVariable
 class Distance:
     """A Distance represents a known distance"""
 
+
     def __init__(self, a, b):
         """Create a new Distance
         
@@ -15,7 +16,7 @@ class Distance:
             b - point variable
         """
         self.vars = (a,b)
-
+    
     def __str__(self):
         return "dist("\
             +str(self.vars[0])+","\
@@ -63,6 +64,12 @@ class Angle:
 
 class Cluster(MultiVariable):
     """A set of points, satisfying some constaint"""
+    
+    staticcounter = 0
+
+    def __init__(self):
+        Cluster.staticcounter += 1
+        self.creationtime = Cluster.staticcounter
 
     def intersection(self, other):
         shared = Set(self.vars).intersection(other.vars)
@@ -125,7 +132,8 @@ class Rigid(Cluster):
         
            keyword args:
             vars - list of variables 
-        """ 
+        """
+        Cluster.__init__(self)
         self.vars = ImmutableSet(vars)
         self.overconstrained = False
 
@@ -153,6 +161,7 @@ class Hedgehog(Cluster):
             cvar - center variable 
             xvars - list of variables
         """ 
+        Cluster.__init__(self)
         self.cvar = cvar
         if len(xvars) < 2:
             raise StandardError, "hedgehog must have at least three variables"
@@ -182,6 +191,7 @@ class Balloon(Cluster):
            keyword args:
             vars - collection of PointVar's
         """
+        Cluster.__init__(self)
         if len(variables) < 3:
             raise StandardError, "balloon must have at least three variables"
         self.vars = ImmutableSet(variables)
diff --git a/geosolver/geometric.py b/geosolver/geometric.py
index 59a9610..c5ba04e 100644
--- a/geosolver/geometric.py
+++ b/geosolver/geometric.py
@@ -328,9 +328,9 @@ class GeometricSolver (Listener):
                 map[geocluster].append(drcluster)
             
         for geocluster in geoclusters:
-            # pick drcluster with fewest solutions
+            # pick newest drcluster 
             drclusters = map[geocluster]
-            drcluster = min(drclusters, key=lambda c: len(self.dr.get(drcluster)))
+            drcluster = max(drclusters, key=lambda c: c.creationtime)
             # determine solutions
             solutions = self.dr.get(drcluster)
             underconstrained = False
diff --git a/test/test.py b/test/test.py
index 1ebdd8d..895bd0d 100644
--- a/test/test.py
+++ b/test/test.py
@@ -739,11 +739,11 @@ def test(problem, use_prototype=True):
     #diag_select(".*")
     print "problem:"
     print problem
+    print "use_prototype=",use_prototype
     solver = GeometricSolver(problem, use_prototype)
-    #solver.set_prototype_selection(use_prototype)
-    #print "drplan:"
-    #print solver.dr
-    #print "number of top-level rigids:",len(solver.dr.top_level())
+    print "drplan:"
+    print solver.dr
+    print "top-level rigids:",solver.dr.top_level()
     result = solver.get_result()
     print "result:"
     print result
@@ -751,6 +751,7 @@ def test(problem, use_prototype=True):
     check = True
     if len(result.solutions) == 0:
         check = False
+    diag_select("GeometricProblem.verify")
     for sol in result.solutions:
         print "solution:",sol
         check = check and problem.verify(sol)
@@ -758,7 +759,6 @@ def test(problem, use_prototype=True):
         print "all solutions valid"
     else:
         print "INVALID"
-
   
 # ----- what to test today -------