Two parts to this,
GIT-related stuff, patches and commits:
Regardless of the correctness of the code, I can't pull in a commit like this in OOFEM main repo. At a quick glance,
1. You are committing a bunch of additional files (Debug?)
2. Non-english comments
3. Unrelated changes (drill stiffness)
So, i couldn't even cherry-pick this particular commit, much less do a pull. Even if you add additional cleanup-commits, I would be to reluctant to pull in all the history.
----------
Now, the computations themselves.
Anything is not better than nothing here. If there is no true, mathematically valid, comparison to make, then we shouldn't allow for such a comparison to be made.
Normal shell elements have some strain-like measure typically called the curvature which is defined from the nodal unknowns and the shell kinematics. For simpler elements, this might be a real, easy to interpret, and comparable quantity. For more complicated kinematics, this might lack a physical meaning. So curvatures of different elements should not be compared.
The shell moment is just whatever happens to be work-conjugate with this curvature. It does not necessarily correlate with equilibrium.
We have 2 paths here. Either we make use of the shell kinematics themselves to obtain sensible definitions of curvature and moment.
Or we could just decide the define the shell moment from equilibrium (and make it clear that this is what is exported).
The curvature is trickier, but we could define is as a linear strain variation over the thickness, and then just fit it to the strain values (since there are only 2 strains, it would produce a perfect fit). As opposed to the moment (which is always sensible defined), this definition would lack meaning if there was more GPs that had a nonlinear relationship, which is why this is less obvious.
I would prefer to see the shells kinematics used to compute this. Taking the GP's natural coordinates and compute what would have been the "curvature" and "moment" according to this particular shell-kinematics. It would basically be like computing a modified B-matrix.
This is difficult to do, since all of mitc4.C is a bit of a hardcoded mess.
If we go for the route of just defining a moment arbitrarily, the sensible choice is of course the integral over the thickness:
M_ij = integral z * sigma_ij dz
Numeric integration over the thickness for the 2 GPs. So, for the ShellMoment, your code is definitely wrong.
Numerical derivative for the curvature.. well, I guess I won't be much better than that. But it's important to know what this is just an a post-processing, and it has nothing to do with the assumed shell kinematics.