Now I understand what was meant by multiple load cases, sorry for the previous confusion.
The mindlin element isn't working with the linearstatic solver, so nonlinear solver is the only option.
And, luckily, multiple elastic load cases are just as well supported in the nonlinear solver.

For the nonlinear solver, I realized that refloadmode isn't set to incremental as default, so you'll also need

``NonLinearStatic nsteps 1 controlmode 1 rtolf 1.e-8 nmodules 1 refloadmode 1``

(I think it should be default, and in fact, with no option to switch it)

This setup isn't smart enough to keep the old tangent for each new load though, but we can tell nlinearstatic to keep the initial secant all steps by setting "stiffmode 3"

``NonLinearStatic nsteps 1 controlmode 1 rtolf 1.e-8 nmodules 1 refloadmode 1 stiffmode 3``

which will keep the same tangent for all steps, along with factorizations.

Mikael,
I still think there is some misunderstanding on the multiple load cases on your side.
The analysis with multiple load cases means, that the simultaneous analyses are made fot the same structure with multiple, independent load cases.
I am pretty sure, that nonlinear analysis does not support this, as this would imply keeping and storing sepate history for each load case, for example.
This multiple load cases are usually suported only by linear analysis, as one can reuse the stifness and there is no need to keep individual histories, one just needs to assemble different load vectors.

In what you have proposed, the solutions in individual steps are not independent, the results are accumulated, as the nonlinear solver incrementally accumulates the solution, no matter how the load is specified.

Borek,

The nonlinear solver will accumulate the loads back to zero after the function peak.
In the file i posted, i do this

``````ConstantPressureLoad 2 loadTimeFunction 1 loadtype 6 ndofs 1 Components 1 100.
ConstantFunction 1 f(t) 1.0
PeakFunction 2 t 1.0 f(t) 1.0
PeakFunction 3 t 2.0 f(t) 1.0``````

which is supposed to produce the total loads: 100, 150, -50.

and with refloadmode = 1 (i.e. the increment == the actual increment) this will produce a behavior like

``````f_1 =       1.0*100. +    0.0*50. + 0.0*(-150.) = 100
f_2 = f_1 + 0.0*100. +    1.0*50. + 0.0*(-150.) = 150
f_3 = f_2 + 0.0*100. + (-1.0)*50. + 1.0*(-150.) = -50``````

stiffmode 3 also retains the tangent matrix over timesteps, so it is functionally identical to the linear solver.

I tested it before posting and everything seems to be correct. The midpoint deflection is
-6.54e-6
-9.82e-6
+3.27e-6
which is the expected values.

Hi Borek, Mikael,
I mange to run my first example which was completely prepared by CAFE (for plate element). Now question is how to calculate element/nodal stresses based on output we have. Is there any option for output so that it gives nodal/elemental stresses in output instead of internal stresses/strains.
Thanks in advance and kind regards, Svemir.

Hi Svemir,

Exporting to VTK, you have the option of "vars" and "cellvars" for all internal states.
Choosing "cellvars" produces one value per cell, and it is computed by a weighted average of the Gauss points.

Choosing "vars" relies on one of the smoothing algorithms available. The simplest one is nodal averaging. Supporting this is optional for all types of elements, and it is not done by this mindlin element yet.

I also think "IST_ShellForceMomentumTensor" is an odd thing to export. It's not a tensor. It's a vector of forces and moments (not momentum). I think it would be far nicer to export forces as a vector, and moments as actual vectors in VTK.
IST_ForceVector
IST_MomentVector
For averaging, we must also consider local coordinat systems. I think CCTPlate (and all other beams and shells) are wrong in the simplistic way they export this data, especially when it is averaged in the nodes (differnet local c.s. == nonsense average)

I can make time for hired work if you wish to have these export features implemented quickly, otherwise you can wait until someone gets around to implement these missing features, or send in a patch yourself.

Regards, Mikael

Hi Mikael,
thank you for your answer. My opinion is that hiring someone who knows code better then we is best choice. Now we are testing complete system because we are preparing application for some projects which will include OOFEM. On that problem we will have money for new development. If project doesn't go we will see can we (BVB CAFE) pay for additional development.

On think in general is not clear to me: how tr_shell element supports all and gives good results while Quad1Mindlin3d you thnik is not promising. They should have same theory in background and Quad must have better results. Can you point to chapter in book you used to create QuadMindlin3d element.
Thanks in advance and kind regards, Svemir.

Kind regards, Svemir.

Hi Svemir,

I'm not an expert on plates. I only know they tend to have tricky locking mechanisms that gives unreliable answers.
The reference book i have doesn't include anything on triangular elements, so I'm not sure how it relates to shear locking.
With reduced integration, the quad element works great, but without any form of stabilization in place they are unrealiable for general usage.
Also, have we checked to see that the CCT element really is that good for thing structures? It doesn't do any reduced integration from what i can see (though perhaps it doesn't exhibit as much shear locking).

I found a digitial copy of the book here
ftp://ftp.demec.ufpr.br/disciplinas/EME … 20733p.pdf

section 15.3  (page 556 of the pdf itself) is where it starts.
Table 15.3-1is useful, as well as figure 15.3-2

Hi Svemir,

I would not be so sceptical about the quality of the element, the reduced integration improves the behaviour for thin plates, but of course, there are better elements for thin plates. The implementation of nodal smoothing procedure is not a complex task, however, at this time I am quite busy with other things. I also have a master student working on several plate elements, but this can take a longer time before being implemented in oofem and checked. So you can wait or try it by yourself, there are many examples already available.
Regarding the CCT element: this is a constant curvature element, which comes with improved lateral defection interpolation (quadratic over the edge) so that the shear is constant. This has the similar effect as reduced integration, but here it is achieved by a smart choice of interpolation functions. So not surprised, when results similar to quad with reduced integration. If interested, I can provide some documentation.

Borek

Hi Borek,
I am not sceptical about those elements at all since we also used reduced integration and different shear coefficients in practical implementation. At the moment the only problem is what sometimes I would like to go bit more faster than I can in real life what means I need a bit more time to understand the code in order to start adding something there (in my local copy). Prior to doing any changes to elements I have to be able to read results properly so that we can compare it with tools that we are using and see if there is any difference at all. From that point of view any documentation you maybe have is welcome. First problem I am faced with is how to calculate Stresses out of results we have (if you can direct me to code which is doing that for existing shell elements). Second problem is related to vtk module. Where I can read what settings are available for vtk output. Maybe I did not read manuals carefully but I still did not find where it is written.
Thanks in advance and kind regards, Svemir.

Borek, Svemir,

I'm surprised you wouldn't be skeptical when using reduced integration with no checks whatsoever if there are spurious modes in the solution.
The selective integration definitely leads to additional zero eigenvalues in the element stiffness matrix.
This is a very well documented problem, and all commercial codes have features implemented to inhibit these spurious modes.
E.g. http://www.dynasupport.com/howtos/element/hourglass
and OOFEM has no such control implemented.
In most cases you can probably get away with it without exciting any of the hourglass modes, but that seems risky without doing any control.

svemir,
Documentation for VTKXML is available here
http://www.oofem.org/resources/doc/oofe … ode33.html

In the input record, you would specify some internal state to export, specified in the header file
src/oofemlib/unknowntype.h

The shell elements only know how to export the shell forces and moments. If necessary, it would be possible to add a new internal state type, such as computing the maximum tensile stress, or something along those lines implement support for that in each shell element.
Currently, there is no shell element that actually does that.

There is a bug in this quad element related to the rather confused definition of IST_ShellForcesMomentumTensor (since it is not really a tensor) so that needs to be fixed before you can get anything useful as nodal values in VTK. The results are still printed in the *.out file however.

I have just added support for nodal value recovery on Quad1MindlinShell3d element, so now you can visualize the moments and forces on this shell element in paraview. When you run the example posted before (svemir01.in) and then in paraview instead of solid color select IST_ShellForceMomentumTensor and partcular component, you can see the results for particular variable.
The element can export strains/curvatures (IST_ShellStrainCurvatureTensor) and forces/moments (IST_ShellForceMomentumTensor)
The component order:
o) IST_ShellStrainCurvatureTensor: eps_x,eps_y,gamma_xy, kappa_x, kappa_y, kappa_xy, gamma_zx, gamma_zy
o) IST_ShellForceMomentumTensor: Nx, Ny, Vxy, Mx, My, Mxy, Vxz, Vyz

Edit: these quantities are in local coordinate system of element, so on planar domains, where the lcs is same for all elements, its ok, while on curved surfaces this will be a problem and one has to implement transformation into global cs.

I have also realized, that each element stores local coordinates of each node. It always difficult to say, whether to prefer speed or memory, but in this case I would prefer to compute local coordinates when needed. Particularly, when large models are considered.

Hi Mikael, Borek
I am testingQuad1MindlinShell3d and it seams that in IST_ShellForceMomentumTensor: Nx, Ny, Vxy, Mx, My, Mxy, Vxz, Vyz values for last Vxz and Vyz are reversed. Can it be related to one of your posts (see below). I was looking for place where it is done but without success. If it is not a problem please letm me Know where it is and I will test it on my copy. Thanks in advance, Svemir.

nabla = [
x  0  0  0  0
0  y  0  0  0
y  x  0  0  0
0  0  0  x  0
0  0  0  0  y
0  0  0  y  x
0  0 -y  0  n
0  0 -x  n  0
]
eps = nabla * [D_1, D_2, D_3, R_1, R_2]

So:
eps = [
eps_11
eps_22
gamma_12
kappa_11
kappa_22
kappa_12 + kappa_21
-eps_32 + R_2
-eps_31 + R_1
]

It is a very real possibility that I made a mistake putting in the components in the B-matrix.
lines: 260-279
(I think you are correct)

Also, when exporting, please note that the smoothing will not take care to respect the elements local coordinate system, so if a quad element is rotated 90 degrees, you will end up with garbage values.
I strongly believe the best way for us would be to add
IST_MomentTensor
IST_ForceTensor
which should contain

``````IST_MomentTensor = [
Mx  Mxy 0
Mxy My  0
0   0   Mz    <---- (drilling-moment)
]``````

and

``````IST_ForceTensor = [
Nx  Vxy Vxz
Vxy Ny  Vyz
Vxz Vyz 0
]``````

both transformed to the global coord.sys.
Then one can use plugins in Paraview to compute eigenvalues and do any kind of post-processing.

If one exports in the element-local-c.s, then it is safe to stick to the local coordinate system (though one should be aware of that).

( Some things could theoretically be smoothed in local c.s, if you think of the local coordinate systems in neighboring elements as "almost" the same. But we have no way of ensuring that. Even if all elements are in a plane, one of them could have its nodes permutated )

Hi Mikael,
is this part that has to be changed?
// shear strains
answer(3 + 3, 2 + 0 + i * 5) = -dns(i, 1);
answer(3 + 3, 2 + 2 + i * 5) = ns(i);
answer(3 + 4, 2 + 0 + i * 5) = -dns(i, 0);
answer(3 + 4, 2 + 1 + i * 5) = ns(i);

Thanks, Svemir.

Yes, that should be it. If I'm not mistaken, the rows should change place;

``````        // shear strains
answer(3 + 3, 2 + 0 + i * 5) = -dns(i, 0);
answer(3 + 3, 2 + 1 + i * 5) = ns(i);
answer(3 + 4, 2 + 0 + i * 5) = -dns(i, 1);
answer(3 + 4, 2 + 2 + i * 5) = ns(i);``````

Mikael,
now it works fine. Please change it in master copy. Thanks, Svemir.

Fixed

Hi Mikael, Borek
I prepared one model with triangle shell element and can not run it. Can you please check input file and tell me what is wrong. Thanks in advance, Svemir.

#*******************************************************************************#
#*****************************BVB CAFE OOFEM EXPORT*****************************#
#****************************    VERSION 01/26/14   ****************************#
#*******************************************************************************#
#OUTPUT FILE RECORD
plate-d.out
#JOB DESCRIPTION
plate-d NonLinear OOFEM
#ANALYSIS TYPE
NonLinearStatic nsteps 1 controlmode 1 rtolv 1E-08 rtold 1E+08 nmodules 1
vtkxml tstep_all domain_all primvars 1 1 vars 1 9 stype 1
#DOMAIN TYPE
domain 3dShell
#OUTPUT MANAGER
OutputManager tstep_all dofman_all element_all
#COMPONENTS SIZE
ndofman 12 nelem 10 ncrosssect 1 nmat 1 nbc 2 nic 0 nltf 2 nset 2
#NODES
node 1 coords 3 0 0.125 0
node 2 coords 3 0 -0.125 0
node 5 coords 3 0.05 0.125 0
node 6 coords 3 0.05 0.075 0
node 7 coords 3 0 0.075 0
node 8 coords 3 0.05 0.025 0
node 9 coords 3 0 0.025 0
node 10 coords 3 0.05 -0.025 0
node 11 coords 3 0 -0.025 0
node 12 coords 3 0.05 -0.075 0
node 13 coords 3 0 -0.075 0
node 14 coords 3 0.05 -0.125 0
#END NODES
#ELEMENTS
#Element not used in oofem 1
#Element not used in oofem 2
#Element not used in oofem 3
#Element not used in oofem 4
#Element not used in oofem 5
#Element not used in oofem 6
tr_shell01 7 nodes 3 2 14 12 boundaryLoads 2 1 1
tr_shell01 8 nodes 3 2 12 13 boundaryLoads 2 1 1
tr_shell01 9 nodes 3 13 12 10 boundaryLoads 2 1 1
tr_shell01 10 nodes 3 13 10 11 boundaryLoads 2 1 1
tr_shell01 11 nodes 3 11 10 8 boundaryLoads 2 1 1
tr_shell01 12 nodes 3 11 8 9 boundaryLoads 2 1 1
tr_shell01 13 nodes 3 9 8 6 boundaryLoads 2 1 1
tr_shell01 14 nodes 3 9 6 7 boundaryLoads 2 1 1
tr_shell01 15 nodes 3 7 6 5 boundaryLoads 2 1 1
tr_shell01 16 nodes 3 7 5 1 boundaryLoads 2 1 1
#PROPERTIES
#PLATES
SimpleCS 1 material 1 thick 0.01 drillstiffness 10000 set 1
#END PROPERTIES
#MATERIALS
IsoLE 1 d 8000 E 209999980051.072 n 0.3 tAlpha 0.0
#END MATERIALS
BoundaryCondition  2 loadTimeFunction 2 values 6 0.0 0.0 0.0 0.0 0.0 0.0 dofs 6 1 2 3 4 5 6 set 2
ConstantFunction 1 f(t) 1.0
ConstantFunction 2 f(t) 1.0
#PROPERTY SETS
Set 1 elements 10 7 8 9 10 11 12 13 14 15 16
#FIXED NODES SET 2
Set 2 nodes 6 1 2 7 9 11 13

Hi Svemir,

thank you, as this has revealed some errors in the code as well (Please update your code to the recent repository version).
Please find attached corrected input file. The main modification is due to the fact, that tr_shell01 element does not support surface loads, but it does support dead weight load, so I changed the input. However, the dead weight load takes into account the material density, so you must take this into account when applying the load (I did not done this, so the example applies  different load value  compared to your original example).

* If you are going to send an input, please use attachments, instead of copying the input in the post, as this makes the thread long and difficult to follow. The attachments can be added when you preview the post or later, by editing it.

Hi Borek, Mikael,
It seams that we will need to have surface load in tr_shell01 element. If I want to add it to tr_shell01 element can I use same logic as it is done in Quad1MindlinShell3d? I was looking in programmers manual without much succes on how to add it. No I would like to try "copy" that part from Quad1MindlinShell3d element. Is it so that both elements have 4 point Gauss Integration? Thanks in advance and kind regards, Svemir.

tr_shell01 is a special element. It only combines two existing elements, and doesn't do much internally.
A surface pressure load would need to be added to "CCTPlate3d",  and a simple wrapper should be added to tr_shell01 and just passes it along to the plate part.

Integration points are not fixed (can be set in input files), but the defaullt is 1 for both elements (since they are linear). But that is just for assembling the internal forces. External forces are free to create whatever integration rule they see fit  (and they do).
See

I have just introduced new IST types for plate and shell elements

``````IST_ShellForceTensor,
IST_ShellStrainTensor,
IST_ShellMomentumTensor,
IST_ShellCurvatureTensor``````

and modified all elements modified accordingly, so that giveIPValue can return these. Modified and unified also the printing methods to separately output forces, strains, curvatures, and moments.

I have also implemented the support for surface load on tr_shell01 element. It is slightly different from Mindlin quad element, as only on one surface (number 1) the load can be applied. The surface load can be specified in local or global coordinate system. Still need extensive testing .....

I think that pressure load would be most suitable here though. Having to specify a traction in the normal direction would complicate the input files a lot.

Also, is there any point to adding NIP 4 ? It isn't needed for the internal forces, and the loads can just as easily set up an integration rule depending on the order of the load instead. E.g. for the linear element and constant load, just 1 integration point is already exact.