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 <alex58spr...@gmail.com> 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 <logar...@googlemail.com>: > > > > 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 <alex58spr...@gmail.com> > 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 > >> >
