Order spheres before points within each coordinate block

In the cases I've tried so far, this leads to substantially smaller
Gröbner bases.

engine-proto/Engine.jl

@@ -118,28 +118,39 @@ equation(rel::AlignsWithBy) = mprod(rel.elements[1].vec, rel.elements[2].vec) -
 # --- constructions ---
 
 mutable struct Construction{T}
-  elements::Set{Element{T}}
+  points::Set{Point{T}}
+  spheres::Set{Sphere{T}}
   relations::Set{Relation{T}}
   
   function Construction{T}(; elements = Set{Element{T}}(), relations = Set{Relation{T}}()) where T
     allelements = union(elements, (rel.elements for rel in relations)...)
-    new{T}(allelements, relations)
+    new{T}(
+      filter(elt -> isa(elt, Point), allelements),
+      filter(elt -> isa(elt, Sphere), allelements),
+      relations
+    )
   end
 end
 
-function Base.push!(ctx::Construction{T}, elem::Element{T}) where T
+function Base.push!(ctx::Construction{T}, elem::Point{T}) where T
-  push!(ctx.elements, elem)
+  push!(ctx.points, elem)
+end
+
+function Base.push!(ctx::Construction{T}, elem::Sphere{T}) where T
+  push!(ctx.spheres, elem)
 end
 
 function Base.push!(ctx::Construction{T}, rel::Relation{T}) where T
   push!(ctx.relations, rel)
-  union!(ctx.elements, rel.elements)
+  for elt in rel.elements
+    push!(ctx, elt)
+  end
 end
 
 function realize(ctx::Construction{T}) where T
   # collect coordinate names
   coordnamelist = Symbol[]
-  elemenum = enumerate(ctx.elements)
+  elemenum = enumerate(Iterators.flatten((ctx.spheres, ctx.points)))
   for coordindex in 1:5
     for indexed_elem in elemenum
       pushcoordname!(coordnamelist, indexed_elem, coordindex)
@@ -167,7 +178,7 @@ function realize(ctx::Construction{T}) where T
   # turn relations into equations
   eqns = vcat(
     equation.(ctx.relations),
-    [elem.rel for elem in ctx.elements if !isnothing(elem.rel)]
+    [elem.rel for (_, elem) in elemenum if !isnothing(elem.rel)]
   )
   Generic.Ideal(coordring, eqns)
 end