For the transportelement, computeBoundaryLoadVector  and  computeBoundaryEdgeLoadVector  are not at all the same though.
1. Different integration rule (edge vs surface)
2. Different "detJ"
3. Different "gcoords"
3. Different "dA/dL"  (at least, i think they should be different, this could be a bug, I don't think the edge load should include the thickness, it might have gone unnoticed since the default is to just return 1.0 here for transport problems, as they don't support CS thickness).

There are of course a bunch of similar code here, but it was difficult to reuse the code, since about every other block of code was different.

What I would suggest, is that we would move more things over to the load.
For example, all this stuff:

        if ( load->giveFormulationType() == Load :: FT_Entity ) {
            load->computeValueAt(val, tStep, lcoords, mode);
        } else {
            load->computeValueAt(val, tStep, gcoords, mode);
        }

        if ( load->giveType() == TransmissionBC ) {
            val.negated();
        } else if ( load->giveType() == ConvectionBC ) {
            field.beProductOf(N, unknowns);
            val.subtract(field);
            val.times( -1.0 * load->giveProperty('a', tStep) );
        } else if ( load->giveType() == RadiationBC ) {
            field.beProductOf(N, unknowns);
            val.subtract(field);
            val.times( -1.0 * getRadiativeHeatTranferCoef(load, tStep) );
        }

should be possible to reduce down to

        load->computeLoad(val, tStep, mode, lcoords, gcoords, unknowns);

Right now, all the code is doing is asking the load for all of these components, and putting the pieces together itself.
This would simplify computeBoundaryLoadVector and computeBoundaryEdgeLoadVector to the point where they almost have no code in common.
It would also allow us to make, for example, nonlinear convective b.c. (though I don't think that's a particularly high priority)
But! there could always be corner cases which I haven't thought about, so trying to separate the load and element like this might not be that simple.

-----------------------------------------------------------------------

I would define the loads as
Nodal loads, N
Edge load, N/m
Surface load, N/m^2
Body load,  N/m^3
regardless of the integration domain.

So, as far as I can think of, three different cases for integrating these types of loads:
Example 1, beam/bar element
A beam element would integrate body loads, and edge loads, over the same edge, but the integration would be different in a single regard: "dx" or "dV" when integrating.
So, this would seem very redundant, since there is so little change, and one could think of reusing the code a lot here.

Example 2, any standard solid element
Say a tet-element, here, edge loads and surface loads are very different. Little to no code-reuse (more than what we already do).

Example 3, 2D/axi-symm elements
Here, one could apply an edge load, which would be on a single node (that goes around the real geometry according to the axi-symmetry.
And the surface load would be defined on one "edge" of the integration domain (and the thickness is analytical the 2*pi*r or from the CS.)

So, right now, we don't have any "EdgeLoad" class type, so it's all the general "boundary", but we need to differentiate if the load is N/m or N/m^2.
That's why I don't think we can merge the two sets into one list.

For some elements, the code might look very similar when integrating loads, but trying to abstract away those differences always seemed like way more hassle than what it is worth. In most cases, we can support many elements with base class implementations anyway, so there are limited implementations to worry about.

But there is definitely code duplication here, I made a new version for assembling loads through sets, and I'll go through the reasons below.
Sorry for another long reply, but there was a lot to go through here.

First off, "boundary" was an attempt from me to keep things general when considering 2D vs 3D elements, but clearly this was failure.
It is probably a good idea to rename them.

The reason why I avoided using the old code was because I didn't want to risk breaking anything. The base of these "new" methods are in the base class Element, and I needed to make this work with Transport and Structural element at the same time.
The plan was always that these should be merged into the same code path eventually, but I needed to be able to test them side-by-side.
I think the transport elements where mixing up what constitutes an edge and a surface, since they don't actually include the thickness of the cross-section anyway. I could be wrong about the intentions there, since if the thickness is 1.0 it will be identical.

These are the three assemble methods that I developed for the input format using sets. (i.e. elements do not contain the record of loads, instead the loads contain a set of elements):

A. computeLoadVector  -  [N/m³] integrates given body load. I.e. deadweight on a beam, that would still be a body load, even if it varies over an edge.
Output for the lspace element is a 24 length array.

B. computeBoundaryLoadVector  -  [N/m²]  integrates given surface load. Loads are defined per area  N/m^2. 2D elements multiple this by the, possibly varying, thickness.
Output for the lspace element is a 12 length array.

C. computeEdgeLoadVector  -  [N/m] integrates given edge load. Typically only for beams, bars, or shells/plates. E.g. a plane strain 2D element would not use edge loads. I attached a little diagram trying to show what I mean here. Surface loads as "t" and edge loads as "q", applying to plane stress and axi-symm elements.
Output for the lspace element is a 6 length array.

The integration for each of these routines is different (different integration rules, calls the interpolator in different ways, different size output etc.)
Trying to reuse code in between these is probably not worth it, they are only superficially similar.

---------------------------------------------

Then there are the "old" methods for load integration (computeEdgeLoadVectorAt and computeSurfaceLoadVectorAt).
These are primarily different in 2 aspects:

1. The output for lspace is 24 length arrays (padded out with zeros) for all load types, which means they have to add an additional mapping to the internal edge/surface and up-assemble the results (i.e. giveEdgeDofMapping)
I found it unnecessary, so I did my routines for only the surface/edge nodes.

2. When they are assembled, they are expected to be in the element l.c.s.  (because the assembler code will transform all the vectors to global c.s.)
In most cases I could think of, it is easier (or at least just as easy) to just compute the load in global c.s.
Since there was no code existing for element lcs -> global c.s. for *only* edge/surface nodes, I opted to just leave it as global c.s. and let the element deal with any transform that's needed.
We could change this, either way, there are only a few elements with local c.s., so it wouldn't make a huge difference.

There is also the differences on what constitutes a "surface" or "edge" load on 2D elements. I think it's best to define let the units define the behavior here. bulk = N/m³,  surf = N/m², edge = N/m.

In the end, I imagine we would have the following routines

    virtual void computeBodyLoad(FloatArray, BodyLoad, CharType, ValueModeType, TimeStep);
    virtual void computeSurfaceLoad(FloatArray, SurfaceLoad, int, CharType, ValueModeType, TimeStep);
    virtual void computeEdgeLoad(FloatArray, EdgeLoad, int, CharType, ValueModeType, TimeStep);

    // And for convection/springbed type b.c.s 
    virtual void computeSurfaceLoadTangent(FloatMatrix, SurfaceLoad, int, TimeStep);
    virtual void computeEdgeLoadTagent(FloatMatrix, EdgeLoad, int, TimeStep);

But to have a common assembly code, we must decide on how the output is defined:
1. Should we pad them out with zeros?
2. Should they be in global or element local c.s.?
There is clear best choice here. Merits and downsides to either way you choose.

I think we are both in agreement on what each type of load should do, the rest is just really just a discussion on what each type of element should support.

Regarding edge loads on non-3D elements, I would still argue the case that it should support surface loads (as well as the edge loads).
Short version:
1. Axi-symm element: Must support surface loads over the "edge" because otherwise you 2*pi*r?
2. Plane stress/strain: Why not support both in case the thickness varies?
3. Shells and membranes: Here we need the out-of-plane surface(s), and the mid-edges.

I've attached another scan of what i think a sensible edge/surface numbering would be for each of these element types.

Long version (sorry):
Axi-symm elements
For axi-symmetric elements, there is really only 1 reasonable option: surface loads. Surely, any reasonable user input input would be N/m² and have it  computed in the code as:

integral  2*pi*r(s)  * t(s) ds

There really isn't even an option here. When would you ever want to supply q(s) = 2*pi*r(s)*t(s)  directly? It would be very hard  to input the very simple case of a constant external pressure if you had to add manually the varying radius as "q(s) = 2*pi*r(s) * pressure", to the point where almost every user would just give up.

Plane stress/strain elements
Usually the thickness is constant, so it's trivial to compute either traction or line load.  But what about varying thickness?
But we could just support both of these options of input trivially. Edge load? Thickness already included. Surface load? Multiply by thickness. I'm not saying that plane elements should include out of plane effects, so, the traction (like all other inputs) may still only vary in the x-y plane.
Thanks to OOP, there would only be one implementation anyway, so why wouldn't we allow it?

E.g. an internal pressure of a boiler cross-section (using plane strain/stress), the only natural choice would be to directly supply the traction "constantpressureload" and have the code compute

integral thickness * pressure * normal ds

Having the option to input the real pressure, as defined by your problem, into the input file directly is surely preferable.
I can't think of an any real scenarios where i would prefer to supply the edge load "q = p*h" as the input.

Shells

We already have analysis running in OOFEM where we have defined the 2 surfaces as the inner and outer surfaces. This is needed to apply tractions to cause delamination when doing XFEM, as the loads are not applied on the midsurface (but rather on each side of a layered composite, pulling the outer layers apart).

Here, i think the "outside" and "inside" are two well defined surfaces that we definitely need. Simpler shell elements may treat this as the midsurface (ignoring the slight thickness offset), so they just define the sign of the normal.

Supporting edgeloads on the midedges around the element is also sensible (like normal).