Re: [Getfem-users] unexpected non-zero values in a stiffness matrix K: a bar (or truss) problem

[Konstantinos Poulios]

Dear Alex,

Your solution "A * E * Grad_u(1,1) * Grad_Test_u(1,1)" will work for horizontal beams. For general beams though, you need to use a unit tangent vector for each element, lets call it T, then you will have an _expression_

"A * E * (Grad_u*T).(Grad_Test_u*T)"

T needs to be calculated beforehand and added to your model as data, either on a MeshFem or on a MeshImd. The same is also the case for AE which can be different from element to element.

I do not have a good suggestion on how to calculate T (in the reference configuration) easily. At the "faces" of the element, which for a beam element are its end points, the weak form language gives you access to T through the "Normal" keyword, but you have to transfer this information to the interior of the element somehow. I can think of some ways, but none of them is very straightforward.

Best regards

Kostas

 

On Tue, Sep 24, 2019 at 10:17 PM Alex Spring <address@hidden> wrote:

Dear Kostas,

Thank you for your quick response. I appreciate your help.
After having your comment, I worked on deriving the weak form language
_expression_ for a bar problem.

The _expression_:
"A * E * Grad_u(1,1) * Grad_Test_u(1,1)"
was derived for a bar with cross-sectional area A and Young's modulus E.

Attached .py Python3 script validates the _expression_.

For those who interested, I wrote some documents on the weak form and
its _expression_ in GetFEM++.
Attached .pdf documents describe a bar and a general problem of linear
elasticity.

Attached files are:
- bar_expressions_compared.py: for comparing stiffness matrices come
from different expressions.
- CheatSheetForGradDiv.pdf: notes for gradient and divergence of a
scalar, vector and matrix.
- LinearElasticity.pdf: notes for a general problem of linear elasticity.
- BarProblem.pdf: notes for the weak form of a bar problem.

*Please let me know if anyone finds errors, corrections or has recommendations.
*I hope the files will be of some help.

Best regards,
Alex Spring

2019.9.19(Thu.) 22:18 Konstantinos Poulios <address@hidden>:
>
> Dear Alex,
>
> Thanks for your interest in the framework. I guess it is all about the use of
>
> asm_linear_elasticity
>
> which was probably never meant to be used in 1D. The brick that you have used is equivalent to the generic weak form language _expression_
>
> (lambda*Trace(Grad_u)*Id(qdim(u)) + mu*(Grad_u+Grad_u')):Grad_Test_u
>
> So basically you need to write the corresponding _expression_ for a bar and use it with the
>
> asm_generic
>
> function instead of asm_linear_elasticity.
>
> Best regards
> Kostas
>
>
> On Thu, Sep 19, 2019 at 2:17 PM Alex Spring <address@hidden> wrote:
>>
>> Dear All,
>>
>> I am handling a bar in GetFEM++ framework. This post is questioning about unexpected non-zero values in a stiffness matrix K. I believe that thinking about it helps deeper understanding of the framework.
>>
>> The problem considered is:
>> """
>> Bar:
>>     0----1
>>     (point: 0, 1; convex: 0 (0-1), segment or 1d-simplex))
>>     (2 dof (u, v) on a single point)
>> Coordinate system:
>>     ^ y, v, Fy
>>     |
>>     ----> x, u, Fx
>> Expected for the convex 0 (0-1):
>>              K          U    =    F
>>     [ 1,  0, -1,  0]  [u0]      [Fx0]
>>     [ 0,  0,  0,  0]  [v0]      [Fy0]
>>     [-1,  0,  1,  0]  [u1]      [Fx1]
>>     [ 0,  0,  0,  0]  [v1]      [Fy1]
>>     (K: stiffness matrix, U: displacement vector, F: force vector)
>>     (A bar should have zero stiffness in the transverse direction (i.e. y- or v-direction here).)
>> """
>>
>> I believe that we can handle a bar using FEM_PK(1, 1) ((1, 1): first 1 for 1d-simplex (i.e. segment), second 1 for 2-nodes on a segment).
>>
>> However, obtained K has unexpected non-zero values:
>> """
>> Obtained for the convex 0 (0-1) :
>>                       K                  U    =    F
>>     [      1,      0,     -1,      0]  [u0]      [Fx0]
>>     [      0,    0.5,      0,   -0.5]  [v0]      [Fy0]
>>     [     -1,      0,      1,      0]  [u1]      [Fx1]
>>     [      0,   -0.5,      0,    0.5]  [v1]      [Fy1]
>>     (K22, K24, K42, K44 are expected to be 0.)
>> """
>>
>> Can we know why K has 0.5 values like this? This behavior is unexpected for me (or possibly for people with structural background), but I guess that this behavior is logically expectable and explainable. I searched getfem-users-archives and looked into FEM_PK, GT_PK and some implementations but could not find the reason. The reason would be about the same interpolation of u and v if |K22|, |K24|, |K42|, |K44| were 1, but actually they are 0.5.
>>
>> Could you please give your comments? Any ideas are welcome. It really helps understanding the framework.
>>
>> Attached is the script for python3. The script explicitly obtains stiffness matrix K using gf.asm_linear_elasticity(...). The other part is similar to Roman's one:
>> [Getfem-users] Solving truss structure in GetFEM++ (via Python)
>> https://lists.nongnu.org/archive/html/getfem-users/2011-03/msg00009.html
>>
>> *I like GetFEM++ way: separative design (GeoTrans, Mesh, Fem, Im, ...); high and low layers (Model, asm_...);  point and convex; and so on.and so on.
>>
>> Best regards,
>> Alex Spring
>>