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

[Yves Renard]

Dear Kostas and Alex,

In fact for a one dimensional domain, even a curved one, Grad_u is a tangent to 
the real element, so that A * E * Grad_u.Grad_Test_u should be correct.

Best regards,

Yves

----- Mail original -----
De: "getfem-users" <address@hidden>
À: "Alex Spring" <address@hidden>
Cc: "getfem-users" <address@hidden>
Envoyé: Vendredi 27 Septembre 2019 20:40:28
Objet: Re: [Getfem-users] unexpected non-zero values in a stiffness matrix K: a 
bar (or truss) problem

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
> >>
>