[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Getfem-users] unexpected non-zero values in a stiffness matrix K: a
From: |
Yves Renard |
Subject: |
Re: [Getfem-users] unexpected non-zero values in a stiffness matrix K: a bar (or truss) problem |
Date: |
Mon, 30 Sep 2019 13:27:43 +0200 (CEST) |
Dear Andriy,
Yes, I thought to a scalar function, sorry. The gradient of a scalar function
is tangent to the real element (it is element_K multiplied by the gradient on
the reference element). But for vector fields, in fact the gradient of each
components is tangent to the real element. So if you want to select the
tangential displacement, you indeed need a tangent vector as you did using
"Normalized(element_K)", yes. May be a small simplification is possible
considering that the gradient of the tangent direction is already tangent to
the element (i.e. using Normalized(element_K) only once) but your formulation
is correct and avoid the potential sign problem.
Best regards,
Yves
----- Mail original -----
De: "andriy andreykiv" <address@hidden>
À: "yves renard" <address@hidden>
Cc: "logari81" <address@hidden>, "Alex Spring" <address@hidden>, "getfem-users"
<address@hidden>
Envoyé: Lundi 30 Septembre 2019 10:34:36
Objet: Re: [Getfem-users] unexpected non-zero values in a stiffness matrix K: a
bar (or truss) problem
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
> > >>
> >
>
>