There are two major changes: a cleanup w.r.t. to the DOWB-matrix
stuff: the DOWB-matrices have eaten the "ordinary" matrices (they
contained the scalar case anyway). This shrinks the needed number of
data-structures and function proto-types (e.g. there is now only one
DOF_MATRIX and one OPERATOR_INFO, one oem_solve_d(), but no longer an
oem_solve_dowb() etc.). This simplifies stuff a lot. Hopefully the
added overhead for the ordinary scalar case is neglectible.

Second main-line of changes: the boundary _interpretation_ (Neumann,
Dirichlet etc.) is now clearly separated from the geometric boundary
classification. Boundary segments are assigned a tag as usual, but
that tag has to be interpreted by the application. Typically the
application now passes a bit-mask to any functions which needs to deal
with boundary conditions, one bit for each boundary type which has to
be taken into account (e.g.: if Dirichlet boundary conditions are
desired for a the boundary segments of type 1 and 255, then the
application would pass a bitmask BNDRY_MASK to

dirichlet_bound(..., BNDRY_MASK, ...);

where BNDRY_MASK has bit 1 and 255 set (and bit 0 to denote that the
mask belongs to a boundary at all).