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

[Andriy Andreykiv]

Dear Yves,

I'm really wondering if Grad_u is really tangent to the real element. Imagine that you move only one node lateral to the beam.

The gradient in global coordinates won't be aligned with the beam, will it?

What we do for Euler beams, for example, is project the hessian of the displacement in the tangent direction using element_K, which is the gradient

of the geometric transformation. I'm sure you can do the same with the gradient.

Here is part of our code (give correct result, of course)

ga_workspace workspace;

auto set = [&](const auto &s1, const auto &s2) {workspace.add_macro(s1, s2);};

//beam director vector
set("d", "Normalized(element_K)");

//Hessian of U (curvatures) in d direction
 //(this is still a 3D vector of curvatures projected on global directions)
set("HessUd", "(HessU . d) . d");
//Hessian of the test functions in d direction
set("HessUTestUd", "(HessTestU . d) . d");

Best regards,

                       Andriy

On Sat, 28 Sep 2019 at 13:28, Yves Renard <address@hidden> wrote:

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