Thanks Mikael. I did some tests and unfortunately there is still something wrong:

1. at node 1 there should be Vz=L*p=25e6 N (actually there is half of that), where L=2500mm and p=-10000 N/mm, and Myy=p*(L^2)/2=-3.125e10 Nmm (actually, Myy=0 and Mzz is not null at the free end). See attached diagrams.

2. is it possible to have the frame forces inside the beam (eg., in a "station" in the middle) or OOFEM would give only the forces at nodes? I'm looking how to modify the code; the alternative is mesh each beam in several parts.

Giovanni

2. the diagrams are obtained with NextFEM Designer that reads the OUT file, not from VTK file. I'll try with the Gauss points: if I understood well, the GP output is only available in VTKXML format? or it can be obtained from the usual output manager?

I see that beam3d has 3 GPs, and they'd be enough for my purposes.
Where do I get the forces? I suppose there would be something  like:

    for (GaussPoint *gp : *this->giveDefaultIntegrationRulePtr()) {
        fprintf(File, "\n  Gauss points forces  at  ", gp->giveNaturalCoordinates(), " ");
        for (auto &val : gp->giveIntegrationRule. ... ) {
            fprintf(File, " %.4e", val);
        }

in Beam3d :: printOutputAt(FILE *File, TimeStep *tStep), right?

Thanks again.
Now I see, we have only stresses and strains. It could be convenient to have the "intermediate" frame forces at GPs. There is any way to easily calculate them? Something like giveForceVectorAt(GaussPoint), I found that some code is available on giveInternalForcesVector, that for now gives only the SUM for all the GPs.

I know there are still the other errors related to the load vector, but actually I get the shear value at each GP (-3.7510e+007) different from the one at the end:

beam element 1 :
  local displacements  0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 0.0000e+000 -5.2083e+002 0.0000e+000 0.0000e+000 0.0000e+000 -2.7778e-001
  local end forces     0.0000e+000 1.2500e+007 0.0000e+000 0.0000e+000 0.0000e+000 2.6042e+010 0.0000e+000 -1.2500e+007 0.0000e+000 0.0000e+000 0.0000e+000 5.2083e+009
  GP  1.1  :  strains  0.0000e+000 -2.1672e-020 0.0000e+000 0.0000e+000 0.0000e+000 -2.4021e-004
              stresses 0.0000e+000 -3.7510e+007 0.0000e+000 0.0000e+000 0.0000e+000 -2.2520e+010
  GP  1.2  :  strains  0.0000e+000 -2.1672e-020 0.0000e+000 0.0000e+000 0.0000e+000 -1.1111e-004
              stresses 0.0000e+000 -3.7510e+007 0.0000e+000 0.0000e+000 0.0000e+000 -1.0417e+010
  GP  1.3  :  strains  0.0000e+000 -2.1672e-020 0.0000e+000 0.0000e+000 0.0000e+000 1.7988e-005
              stresses 0.0000e+000 -3.7510e+007 0.0000e+000 0.0000e+000 0.0000e+000 1.6864e+009

while it must have a value that (in the correct case) drops linearly to zero, or (in the actual case, with the error in the distributed load) must have the same value as in both end (1.2500e+007). Am I wrong?

Ok, so the end forces of the beam isn't giving such useful values here, due to the external forces being applied with a set.
In this new way to apply b.c.s, the elements doesn't keep track of their list of loads (which is a good thing, we map loads -> set -> elements, as opposed to the elements->loads)

This means that instead of outputting the end forces, the beam just exports the internal forces.
The element reports "1.2500e+007" in the fixed end, and "-1.2500e+007" in the free end. Add the contribution (2 * 1.2500e+007) of the load, and you'd end up with the expeted 2.5e+7 and 0 respectively.

I can't think of a good way to deal with this. Beams are very much the exception here, as no other element does this AFAIK (and only 1D elements can do this, so beams and bars only)

In short, I think the load is correctly computed, and it all seems to be correct (only it's the internal forces, and not the node forces you get in the output file).

A few points here:
1.
edubeam is just specialized for beam elements, and we have to try to keep OOFEM more general. Adding very specialized conditions for beam elements complicates other things. For example, the "local end forces" output. This requires the old way of applying b.c.s, because the element needs to know what loads are applied. The new way of applying loads (using sets) has huge advantages when it comes to things like adaptivity etc., but it means that the beam no longer keeps tracks of loads, and instead of getting "local end forces = internal forces - external load", you only get "local end forces = internal forces - 0". This is not intentional or desirable, but I don't see a nice way to fix that, because I like the advantages of using sets. I'd rather just output the displacements for the beam.
Nodal loads are not distributed over the element, so these are not part of the external loads on the element.
All this behavior is completely expected.

2.
Why do we want the end forces? Well the Euler-Bernoulli type beam kinematics doesn't satisfy the equilibrium, so we "cheat" and post-process the shear stresses from solving the simple 1D equilibrium.
This is because we don't get the relevant values from the GPs (which are the values that come from the beam kinematics, not equilibrium). They aren't really incorrect.
Now, this post-processing is only possible for beams, and only beams. For a shell, you'd have no option.

3.
Integrating loads in 1D is very simple, so a post-processing tool could technically still use the internal forces for a beam, and just add the external forces. It's more work for the post-processer, but doable. I would bet this is basically how edubeam does it.

So, I've stopped and thought through these issues a bit more carefully now.

1.
I agree with the last part there, there shouldn't be two ways, but the direction I'm in favor off is to remove the  this problem would be to get rid of loads not applied via sets. The end forces are not as useful as one might think (see below)

2.
What I was trying to say was that the post-processing equilibrium is only ever possible in 1D specimens, i.e. beams (and bars).
You can't post-process using equilibrium for shells, because it won't be statically determined.
The forces per unit length can't be obtained (if they were, they would be associated with a an edge, not a node).
This has nothing to do with nonlinear/inelastic solutions, it's just impossible.
For shells/plates, the GP values are the only option.

3.
I'm not sure why you would classify them as linear and nonlinear.
There is no distinguishing needed in the post-processing of beam results. Just take the distributed loads (if any) and integrate those along the element (applying the internal forces at each end).
I would imagine that all FE-software that shows force-moment diagrams for beams are doing precisely this in their visualizers.
If you get the internal forces for the beam, you can just add the distributed load on top of that (gradually integrating over the beam) to obtain the diagrams.
In fact, the old "local end forces" did not give you enough information, because you don't know the shape of the load.
You could also use the moment from the GP-data, and use that to post-process the shear force (which is your only option for shells).

I have made some changes in the B-matrix, so that problem *might* be fixed.
I can't provide an original formulation, because this element dates back to the initial import from Borek (back in 2007), making it at least that old, so at least 3 years before i even started using OOFEM myself.
From what i can tell, it's like an Euler-Bernoulli beam with some shear correction terms (which are there to soften the beam behavior slightly, not to give accurate shear force values).

Yes, axial loads (which is just a bar) is just 1D equilibrium with no kinematic assumptions. The GP values are then just the equilibrium values.
It's just the shear force for beams that you cannot obtain through GPs.
The small deformation beam/shell elements are always just superimposed bar+transverse deformation+torsional axis  & membrane+plate respectively.
The moment should be good in the GPs as well. It's just that one cannot yet the shear forces directly.

Whether one does the post-processing to obtain shear forces inside OOFEM or not remains a question. Doing it with the standard printOutput functionality is quite limiting by its design.
There are several other export modules, and more could be added (like a specialized beam export module).
This is the only (somewhat clean) way to do it.

Hi,

the problem is now solved and fix is submitted to the repo (http://oofem.org/gitweb/?p=oofem.git;a= … 8e8e155e80)
Borek