Here i'm documentating my "cleanup" work. I believe I have some good design changes which will reduce the overall complexity in OOFEM.

What's special about CBS:
The CBS algorithm splits up pressure from velocity, and by doing so it uses EquationID for the split.
In order for that to make any sense, there can only be pressure and velocity dofs in the domain, there can't be any slave dofs or xfem or anything like that (for understandable reasons), otherwise the algorithm doesn't seem to be well defined.
It also introduces AuxMomentumBalance, which isn't actually a separate equation since it doesnt have separate dofs or equation numbers, and it just needed to inform CBS::giveUnknownComponent to grab the auxiliary solution.
Noone has touched the CBS code in quite some time, while the rest of OOFEM has been steadily maturing.
So it's high time to look over this code now since many things can be done better.

What about everything else:
For all other types of problems, we always set up the total system of equations, so we end up passing something like EID_MomentumBalance_ConservationEquation into every single function when we basically mean "just take everything". In fact, since many problems introduce additional lagrange multipliers, slave nodes or xfem, the equation ID isn't even well defined.. So for everything except CBS, EquationID lacks any meaningful definition.
It could be removed with minimal changes, without losing any features. In fact, most code that uses EquationID just leaves a false impression, for example, SUPG, which has

int
SUPG :: giveNumberOfEquations(EquationID id) {
...
    if ( ( id == EID_MomentumBalance ) || ( id == EID_ConservationEquation ) || ( id == EID_MomentumBalance_ConservationEquation ) ) {
        return numberOfEquations;
    } else {
        _error("giveNumberOfEquations: unknown equation id");
    }
...

which gives the impression that it would be possible to assemble just the conversation equations, which it certainly is not.

What should be done:
1. Introducing a different ValueModeType (e.g. VM_Aux)
2. Use a PrimaryField for the auxiliary velocities and have the element pass that to computeVectorOf
3. Introducing V_u_aux ... dofs.
which all seem like reasonable choices.

ConservationEquation and MomentumBalance - split
Instead of using EquationID to separate velocity- and pressure-related matrices and forces, I think a custom UnknownNumberingScheme should be used instead.
E.g. instead of

            this->assemble( mss, stepWhenIcApply, EID_MomentumBalance, MassMatrix,
                           EModelDefaultEquationNumbering(), this->giveDomain(1) );

one could have

            this->assemble( mss, stepWhenIcApply, MassMatrix,
                           VelocityEquationNumbering(), this->giveDomain(1) );

which is a design i think makes much more sense. In fact, EngngModel::giveNumberOfEquation should have a UnkownNumberingScheme parameter.

I've implemented this and it works fine.
AuxMomentumBalance - some odd code
My second step was to rework the code so that it doesn't need EID_AuxMomentumBalance, which was fairly straight forward, except for one part which I couldn't really grasp

void
TR1_2D_CBS :: computeDensityRhsVelocityTerms(FloatArray &answer, TimeStep *tStep)
...
    this->computeVectorOf(EID_MomentumBalance, VM_Total, tStep->givePreviousStep(), u); // THIS WILL INCLUDE DIRICHLET BCs
    this->computeVectorOf(EID_AuxMomentumBalance, VM_Incremental, tStep, ustar); // THIS WILL INCLUDE DIRICHLET BCs
    u.add(theta1, ustar); // This will produce 1.5 times the boundary condition for the prescribed dofs.

so the final velocity u will actually have velocities larger from the actual Dirichlet b.c.
I'm not sure if this was intended, or accidental. You can't really tell from the tiny test examples, so I've left it as is.