Mikael, I forced the V1,V2,V3 and V4 vectors to be the element normals in the code, but I still obtain some errors. Just to know, which will be normal on the example I posted? I gues it will be (1,0,0).

For the poor results: in my case, my aim to check if the results obtained are better than Quad1MindlinShell.

EDIT: I cannot make it work, there are problems related to the stiffness matrix (resulting displacement are all 1.IND#), or it seems so... I remember that you found a bug on other shell elements, see here , as they were tested to work only in XY plane. It seems that this is the same, and there some vectors hard-coded in mitc4.c (for example, e2 in computeBmatrixAt).

Thanks Mikael, I didn't notice your commit in your github. I think I made a similar change in the code, and the same error happens with my copy.
I'm trying to understand what's wrong in computeBmatrixAt, it is not easy for me.

I used so far other shell elements in other solvers made upon MITC4 formulation, but none of them refuse to work with element normal; the behaviour of MITC4shell in OOFEM is not clear to me.

I never looked into why it was singular. It might be as simple as you state Martin.

What do you think about the option of using the z-axis of the nodal LCS for the normal? I think it makes sense, and doing so would aso make it easy to specify the drilling dof for arbitrary orientations. You would need something like that to handle arbitrary geometries.
If you have a T-section, or an intentionally sharp corner, you can just duplicate nodes and use slave-dofs.

Giovanni: You say you get the *same* error? I do not. It passes that part of the code fine, but fails during Newton-iterations as the problem is singular (most likely due to the unrelated issue pointed out by Martin). This would happen to most shell formulations, unless they include a ficticious drilling stiffness.

We could probably add this to the MITC4 element as well.

I can only speculate why you would have to clean the build. We aren't doing anything special in OOFEM w.r.t this. CMake only generates the build scripts (all standard), so I guess the only possible blame is to Visual studio.

WallsB.out
-> single quad model
StaticStructural nsteps 1
domain 3dShell
OutputManager tstep_all dofman_all element_all
ndofman 4 nelem 1 ncrosssect 2 nmat 1 nbc 3 nic 0 nltf 2 nset 2
node 30001 coords 3 3 0.0 0.0 
node 30002 coords 3 3 0.2 0.0 
node 30017 coords 3 3 0.0 0.2 
node 30018 coords 3 3 0.2 0.2 
#Quad1MindlinShell3d 30001 nodes 4 30001 30002 30018 30017  crossSect 2
mitc4shell 30001 nodes 4 30001 30002 30018 30017  crossSect 2
# fictitious section for LumpedMass elements
SimpleCS 1 thick 0 drillstiffness 0 material 1
SimpleCS 2 thick 0.25 drillstiffness 1e11 material 1
IsoLE 1 d 0.0 E 11000 n 0.35 tAlpha 0
BoundaryCondition 1 loadTimeFunction 1 DOFs 6 1 2 3 4 5 6  values 6 0 0 0 0 0 0  set 1
NodalLoad 2 loadTimeFunction 2 Components 3 1 1 1 DOFs 3 1 2 3 set 2
BoundaryCondition 3 loadTimeFunction 1 DOFs 1 4  values 1 0  set 2   # Prescribed the drilling dofs (YZ-plane == R_u dof must be removed)
ConstantFunction 1 f(t) 1.0
ConstantFunction 2 f(t) 1.0
#PiecewiseLinFunction 2 npoints 2 t 2 0.0 1 f(t) 2 0.0 1.0
set 1 nodes 2 30001 30002
set 2 nodes 2 30017 30018

This works fine for me.

I don't disagree about the convenience, and hopefully Edita will look into adding it in the future. But, this type of missing drilling dof control is almost always missing from plate/shell elements.
Look into any textbook, be it a mindlin, kirchoff or MITC4 shell, none of them will mention the drilling DOF. It is just not part of the bending kinematics.
The drilling stiffness "trick" that I added to the quad1mindlin3d shell has nothing to do with mindlin theory.

If you have something like the drill stiffness I added to the mindlin shell, these are concerns:
1. Curved surfaces (or sharp corners, T-sections, etc.) this ficticious drilling stiffness partially adds to the "real" system of plate equations.
2. You run the risk of attaching a different elements, like a rigidarm node, or a beam element, etc. etc. to this "fake" drilling dof, and get answers which are pretty much completely unreliable.
That's not to say that excluding it would be any more *correct* than including it though. Just know that drilling stiffness can affect the solution, and when it does, I wouldn't consider the results reliable.

These equations could also be soved if you had a more forgiving solver. I don't know if it's feasible to create a direct solver that can handle these singularities, but it must be theoretically possible to still solve these types of problems (without drilling stiffness).

I'm simply pointing out that there will pretty much never be an option here where you do not need to know precisely what you are doing (at least not in OOFEM). Every user of OOFEM will need to know what assumptions they are making, understand that they simply get what they ask for, i.e. the elements you ask for is the approximation you get based on FEM.

I'm struggling to integrate shell generalized stresses along a path, let's say along a section cut in a shell wall.
Consider for example this simple model (cantilever with a transveral load at the free end):

three_bis-stV1.out
-> shell model loaded in plane
LinearStatic nsteps 1 nmodules 1
nrm tstep_all domain_all vars 6 9 10 11 104 1 4 stype 1
domain 3dShell
OutputManager tstep_all dofman_all element_all
ndofman 33 nelem 20 ncrosssect 2 nmat 1 nbc 4 nic 0 nltf 2 nset 4
node 1 coords 3 0 0 0 
node 2 coords 3 250 0 0 
node 3 coords 3 500 0 0 
node 4 coords 3 750 0 0 
node 5 coords 3 1000 0 0 
node 6 coords 3 1250 0 0 
node 7 coords 3 1500 0 0 
node 8 coords 3 1750 0 0 
node 9 coords 3 2000 0 0 
node 10 coords 3 2250 0 0 
node 11 coords 3 2500 0 0 
node 12 coords 3 0 0 -250 
node 13 coords 3 250 0 -250 
node 14 coords 3 500 0 -250 
node 15 coords 3 750 0 -250 
node 16 coords 3 1000 0 -250 
node 17 coords 3 1250 0 -250 
node 18 coords 3 1500 0 -250 
node 19 coords 3 1750 0 -250 
node 20 coords 3 2000 0 -250 
node 21 coords 3 2250 0 -250 
node 22 coords 3 2500 0 -250 
node 23 coords 3 0 0 -500 
node 24 coords 3 250 0 -500 
node 25 coords 3 500 0 -500 
node 26 coords 3 750 0 -500 
node 27 coords 3 1000 0 -500 
node 28 coords 3 1250 0 -500 
node 29 coords 3 1500 0 -500 
node 30 coords 3 1750 0 -500 
node 31 coords 3 2000 0 -500 
node 32 coords 3 2250 0 -500 
node 33 coords 3 2500 0 -500 
Quad1MindlinShell3d 1 nodes 4 1 2 13 12  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 2 nodes 4 2 3 14 13  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 3 nodes 4 3 4 15 14  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 4 nodes 4 4 5 16 15  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 5 nodes 4 5 6 17 16  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 6 nodes 4 6 7 18 17  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 7 nodes 4 7 8 19 18  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 8 nodes 4 8 9 20 19  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 9 nodes 4 9 10 21 20  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 10 nodes 4 10 11 22 21  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 11 nodes 4 21 22 33 32  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 12 nodes 4 20 21 32 31  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 13 nodes 4 19 20 31 30  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 14 nodes 4 18 19 30 29  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 15 nodes 4 17 18 29 28  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 16 nodes 4 16 17 28 27  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 17 nodes 4 15 16 27 26  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 18 nodes 4 14 15 26 25  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 19 nodes 4 13 14 25 24  crossSect 2 mat 1 reducedintegration
Quad1MindlinShell3d 20 nodes 4 12 13 24 23  crossSect 2 mat 1 reducedintegration
# fictitious section for LumpedMass elements
SimpleCS 1 thick 0 drillstiffness 0
SimpleCS 2 thick 300 drillstiffness 1e11
IsoLE 1 d 2.5E-09 E 30000 n 0.3 tAlpha 0
BoundaryCondition 1 loadTimeFunction 1 DOFs 6 1 2 3 4 5 6  values 6 0 0 0 0 0 0  set 1
BoundaryCondition 2 loadTimeFunction 1 DOFs 6 1 2 3 4 5 6  values 6 0 0 0 0 0 0  set 2
BoundaryCondition 3 loadTimeFunction 1 DOFs 6 1 2 3 4 5 6  values 6 0 0 0 0 0 0  set 3
NodalLoad 4 loadTimeFunction 2 Components 1 -10000 DOFs 1 3 set 4
ConstantFunction 1 f(t) 1.0
PiecewiseLinFunction 2 npoints 2 t 2 0.0 1 f(t) 2 0.0 1.0
set 1 nodes 1 1
set 2 nodes 1 12
set 3 nodes 1 23
set 4 nodes 1 22

From the .out file, I get the following generalised stresses from elements 1 and 20 (they host the trescribed nodes 1, 12 and 23) in tabular form:

Elem    GP    s12    length L    L * s12
1    1    -4.46E+01    125    -5.57E+03
1    4    8.46E+01    125    1.06E+04
20    1    -2.53E+01    125    -3.16E+03
20    4    6.53E+01    125    8.16E+03
                
            Sum L*s12    1.00E+04  # this is exactly what expected, i.e. the load I set as input

The same results cannot be collected using ZZnodalRecovery. In the nodes 1, 12 and 23 I get:

Elem                
23    sii    length L    L * sii    
s11    -3.53E+02    125    -44124.9982    
s22    -7.70E+01    125    -9621.2503    
s12    5.73E+01    125    7161.2500    
12                
s11    -2.98E-12    250    0.0000    
s22    -6.11E-13    250    0.0000    
s12    4.62E+01    250    11539.9998    
1                
s11    3.53E+02    125    44124.9982    
s22    7.70E+01    125    9621.2503    
s12    5.73E+01    125    7161.2500    
                
        total N    0.0000    
        total V    25862.4997    

# L is the "length of influence" of each node                

# such values are taken from:
    Type       Node                s11        s22           s12
IST_ShellForceTensor    1    3.53E+02    76.9700    5.73E+01
IST_ShellForceTensor    12    -2.98E-12    0.0000    4.62E+01
IST_ShellForceTensor    23    -3.53E+02    -76.9700    5.73E+01  <--------- !!!

What looks strange to me is that I would expect the sign of s12 in node 23 as the opposite than in node 1.
Can you help explaining this? Am I missing something? Sorry for the long post.

EDIT: I'm using NRM_ZienkiewiczZhu recovery model, it is the only one that works with this file.

Mikael, the model should return shear in-plane stresses, their sum cannot be zero despite the rough mesh, right? It seems you're plotting s11.

Also, i don't understand the question.
When doing stress recovery we have, of course, no guarantee that the recovered (and extrapolated) stresses satisfy equilibrium.
The sum? Do you mean the integral? The sum isn't meaningful.

You basically have a 2D, plane-stress, beam, that you load on one end. The plate equation itself does *nothing* in this load case. So you could just as well simplify to use the 2D elements.

So, 2 element, bilinear elements, for computing a beam-type deflection. This is, of course, way way to coarse to capture the true kinematics. Just like you can't get shear stresses in a euler-bernoulli beam, you can't get them from this coarse mesh.
if I'm not mistaken, the "true" variation of s_12 should have a quadratic shape.
So, make a finer mesh. Use VTK to show the results.

It looks like a standard least-square fit to extrapolate the GP values to the nodes. It doesn't really have anything specifically to do with MITC4, and you can find it in some other 4-node elements as well.

It's just fitting this polynomial

answer.at(j) = b.at(1, j) + x1 *b.at(2, j) * x2 * b.at(3, j) * x3 * b.at(4, j);

(though x3 and b.at(4,*) isn't actually needed here, it' just a bilinear interpolation in a plane).

1. The bug:
That seems reasonable. If so, then the bug exists in all other elements that use the same approach (the same code is present in at least the LSpace element).
To be honest, I'm surprised by this bug. Seems like the multiplications should cause huge errors that are easily spotted (but maybe that's only true for very coarse meshes).
I also think the "w" and "x3" parameter must be removed from MITC4, since there only is a midplane here.

2. Local basis vectors:
I think it should only be the computeLocalBaseVectors.

But I don't actually think you should change anything in NodalAveraging...
It's defined in a way that matches the interpolator and integration rule (corner are defined in order: [1,1], [-1,1], [-1,-1], [1,-1]).
If this somehow isn't consistent with the rest of the MITC4 element, then there is another bug.

3. What you get is the integration point values, as they are used for the FE-simulation. By the kinematics assumed by all beams and shell formulations, the shear stresses are (unavoidably) very poorly captured (only weakly captured).

Basically, the only remedy I can think of is
1. Use 3D elements.
2. Some method like  http://www.sciencedirect.com/science/ar … 2315011265 (here the shear stresses are re-solved from equilibrium using the least-square fitted normal stresses)

The "trick" used by Giovanni for beams (see his beam export module that he posted on the forum) is basically just solving each beam analytically. Unfortunately, this only works for 1D models, and cannot be applied to plates.

Also, be aware that the first thing you do with FEM, is to project the loads on the basis functions. All information about the variation across an element is lost. The Gauss points doesn't "know" if the loads come from a uniform pressure, or concentrated nodal forces (they are 100% equivalent in FEM).

The only third alternative is to not look at the shear stresses, because they are completely unreliable.

Hello everyone,
I just wanted to say that I tried by changing the local coordinates for
MITC4Shell :: NodalAveragingRecoveryMI_computeNodalValue
and
FEI2dQuadLin :: giveDerivatives
to match the new local axes and now it works fine.
The "only" thing is that now all other 2d elements and all other FEI2dQuadLin methods may need some refactor as they might be broken because of this.
I'd suggest to modify or overload intepolation functions to get nodal natural coordinates directly from the element or as a function argument.

I agree. But I'd like to highlight the fact that quad1mindlin uses the same interpolation functions with the hardcoded derivatives and stuff but calculates the local axes the way I changed mitc4 (n2-n1, n4-n1). So my guess some results may now be more accurate.
The point is: is it better to use (n4-n3,n2-n3) as local axes for all 2d elements to avoid changing the interpolators?