If they, preferably, explicitly express the "B" matrices, then it will be easy to implement.

I must be clear that none of those changes will ever make it into OOFEM. GeneralBoundaryCondition is a base class for all types of boundary conditions (Neumann/Newton/Dirichlet/???) and should not have an array. Even if you cast to BoundaryCondition, you should still not access the underlying arrays directly. Here, you are not using the boundary condition for anything it was designed for (there are no DOFs here, you ignore the load time function, etc.).
I don't know what you are using the displacement value for so I can't speculate in what I would suggest you do instead, but I'm sure there are alternatives to this approach.

It is quite clear that these are highly specialized customizations that you have done, but to be included in main OOFEM, it must be written in a much more general way. If we don't stick to this it will be difficult to improve OOFEM in the future.

Yes they can be used without Abaqus (it's just a fortran function). There is a very tiny file to include "aba_param_dp.inc" which typically contains 2 lines, one which defines some implicit variables, "implicit real*8(a-h,o-z)" and one that defines whether or not it's a double or single precision machine (is there ever any single precision machines around?)  "parameter (nprecd=2)"

libeam2d is supposed to be a Mindlin beam, so it should be modified. Adding a new formulation would be fine.
beam2d wouldn't be of any use for you either. It's still only valid for small deflections.

For implementing a new element, there are 2 things (minimum) that needs to be supported:

1. Computing the internal forces. In most elements, this is done by "q = integral B^T * stress(B*u) dV". u are the nodal unknowns, and B*u is the strain, evaluated at each integration point. For beams and structural elements, stress and strain refer to the generalized stresses and strains (includes curvatures).

2. Computing the stiffness matrix. In most elements, this is done by "K = integral B^T * D * B dV" by numerical integration. The D-matrix is taken from each integration point.
This is only used for the Newton-iterations (and if you use an alternative solution strategy, you could technically do without K).

So, as long as one can compute the internal forces (q), and the stiffness (K) one can solve the FE-analysis with normal Newton iterations.
One probably also wants to support external loads; "f = integral N^T * t dA", though support for this is optional.

The elements don't need to use the standard B^T * D * B formulation to compute the stiffness matrix. Some elements (like the large deformation shell element "shell7") can't really use this formulation, and if so, then the elements are free to overload the functions to compute the stiffness matrix and do whatever they want.

I would suspect that a large deformation beam element would need to do such a completely different formulation. It is probably a lot harder to implement.

Note on element local c.s.
Also worth adding is the local coordinate systems. Elements can define a local c.s. and compute all the vectors and matrices in this local c.s., which is then transformed into global c.s. automatically. For example, the beam2d and libeam2d elements do this.

For a large deformation structural elements (like the shell7 elements), defining a local c.s. doesn't simplify anything, since it's curved in space. So, typically, the local c.s. is only ever used for small deformation elements.

Do you really need a beam?
Just as a side-note. I'm wondering if you really need beam/shell elements. The beam kinematics and all the assumptions about shear stresses etc. might not be very accurate for forming processes.
All the elements you need for a forming analysis already exists with solid elements.

In my experience, going to textbooks is the only option. Research articles is a dead end, because they typically write things in a speciaized "quiz" way as you describe it

When it comes to OOFEM, work is typically only done while there is a project (often a PhD student). Few research projects involve implementing existing beam models (even though it would be very valuable). Often something new and fancy instead.

Well, in under 30 seconds you could still probably have many tens of thousands of unknowns, so I don't really think that would limit you to much  when it comes to solid elements.  Of couse, it always depends on how thin geometries you have (something like the video you showed earlier would probably be doable in solid elements under 30 seconds).

A rough count, ~30 elements from the edge til the middle,  say, 5 elements in thickness. Solving that problem in 2D (plane stress) should be very accurate.
Using symmetry, and only 1 element in thickness (since nothing really happens in thickness here), just like with your beam examples.
So 31 * 6 nodes, 3 dofs per node, 31 * 6 * 3 = 372 dofs (for a 2D problem), which is nothing. It's hard to tell how many iterations you will need to obtain final converged shape.
If there is some performance sink, I'm sure things could be done to improve that.

I solve system with several hundred thousands of unknowns, and it takes a couple of seconds (partly because I have no lagrange multipliers in my system, so that I can use a iterative (conjugate gradient) solver).

Is it much work for you to duplicate the nodes in thickness for your pre-processor, and add some plane stress elements?
I think this is worth checking out, as beam models for this type of problem can get very very complicated (I browsed through Zienkiewicz book yesterday, and it was over 50 pages long to derive the rod element).

(It's important that you run OOFEM with release mode, and good with a good lapack library and suitable solver, as these options together will speed up OOFEM at least 100 times over compared to a simple debug compilation).

There are plane strain element with nlgeo 1, see for example quad1planestrain.C. However, it depends on what exactly you want/need, since appropriate large strain material model is also needed, and there are only a few large strain models in oofem.

Nlgeo 2 doesn't exist in oofem anymore.

The problem is that the manual is old and need to be updated, sorry for that. Nlgeo 0 means that the element does not support geometrical nonlinearity, i.e. the formulation is based on the "small-strain" stress tensor which is computed form the small strain tensor.

The NLgeo1 is the total lagrangian formulation using the first Piola-Kirchhoff stress and the deformation gradient.  The element needs to implement computeBHmatrix function which computes displacement gradient(in contrast with computeBmatrix which gives its symmetric part), see  structural2delement.C for example of implementation of this function.

If I am right the tutorialmaterial is the vonMises small strain plasticity. If you want to use it in large strain setting it has to be extended to the large strain range using multiplicative split of deformation gradient or some hypoelasto-plastic formulation. Other possibility is to use additive split of logarithmic strain(see lsmastermat.C).

Martin