I've introduced sets, and with these a new method of applying boundary conditions (and as a sideeffect, also initial conditions).
A set is defined by either elements, boundaries of elements, element edges (in 3d) or nodes. Any combination of these is also possible.
Boundary conditions should now be able to simply point towards a set (regardless of Neumann/Dirichlet/Active type)

The old approach for boundary conditions have been split into 3 parts.
1. Neumann boundary conditions. Each individual element kept a list.
2. Dirichlet boundary conditions where applied directly to DOFs.
3. Active boundary conditions kept a list of applied elements.

This new approach using sets unifies the input method for all 3 types (it also follows naturally to add sets to initial conditions as well).
A few things to point out:
1. Internally, DOFs still keep track of their dirichlet BCs and ICs (as they have to). This is just changes the way it's specified in the input file.
2. Domain-type is redundant with this new approach. New code automatically detects all necessary dofs from the elements themselves.
3. RigidArmNodes will need to specify its dofidmask.
4. Whenever one needs to specify doftypes one also needs to specify dofidmask.
5. Dirichlet b.c.s need to specify which dofs they prescribe.

To give you a sample of what input files could look like (old input method left in after # for comparison):

axisymm01.out
Test of linear axisymmetric quadrilateral elements, a thin ring with internal overpressure.
LinearStatic nsteps 1 nmodules 0
OutputManager tstep_all dofman_all element_all
ndofman 6 nelem 2 ncrosssect 1 nmat 1 nbc 2 nic 0 nltf 1 nset 2
Node 1 coords 2 4.0  0.0 #bc 2 0 1
Node 2 coords 2 4.0  0.1 
Node 3 coords 2 4.1  0.0 #bc 2 0 1
Node 4 coords 2 4.1  0.1 
Node 5 coords 2 4.2  0.0 #bc 2 0 1
Node 6 coords 2 4.2  0.1 
L4Axisymm 1 nodes 4 1 3 4 2 NIP 4 crossSect 1 mat 1 #BoundaryLoads 2 2 4
L4Axisymm 2 nodes 4 3 5 6 4 NIP 4 crossSect 1 mat 1
SimpleCS 1
IsoLE 1 d 1.0 E 1000. n 0.0 tAlpha 0.000012
BoundaryCondition 1 loadTimeFunction 1 prescribedvalue 0.0 defaultDofs 1 2 set 1
ConstantEdgeLoad 2 loadTimeFunction 1 components 2 10.0 0.0 loadType 3 ndofs 2 set 2
ConstantFunction 1 f(t) 1.0
Set 1 nodes 3 1 3 5
Set 2 elementBoundaries 2 1 4

Set number 1 could be made out of a list of elementBoundaries. The afflicted dofs on that boundary is automatically detected either way. Switching between a prescribed load and a prescribed displacement is then just a trivial 1 line change.
Sets of elements and nodes can also be conveniently specified by using ranges (input fields "elementRanges" and "nodeRanges") so including all elements for some initial condition is very simple "Set 1 nodeRanges {(1 100)}".

I've implemented backwards compatibility with the old approach, however I think it's a good idea to phase out the old method eventually (it gets a bit messy).
I'm still fixing the last few bugs (making sure all the tests and benchmark problems work) and I'll publish the changes. All old input files should still work.
From here on I'll probably write a python script to convert files from to use the sets instead of the old approach for bcs and ics.
Then when I can verify that everything is working, update the documentation to the new approach.