getfem-users
[Top][All Lists]
Advanced
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Getfem-users] Curvilinear structures in Getfem

From: Jean-François Barthélémy
Subject: Re: [Getfem-users] Curvilinear structures in Getfem
Date: Sat, 2 Nov 2024 18:07:16 +0100
Dear Konstantinos and Yves,

It's been a long time now since I did that and in fact it was just a one-shot try in the framework of a project with students.

From what I can remember, the first point is that it is indeed complicated to define linear local longitudinal displacement of Lagrange type and second order transversal ones of Hermite type due to displacement continuity problems at nodes where the tangent is not continuous (at 90° corner bend for instance at which displacement components change of roles between longitudinal and transverse parts). And furthermore a full Hermite type displacement is not sound either I think due, once again, to problems of discontinuities at corner bends (the continuity of the rotation is expressed by means of gradients with respect to different components of displacements from one side to the other) or even undesired continuities: for instance we shouldn't impose the continuity of the gradient of longitudinal displacement (wrong at junction between two materials).

The only simple sound way I found was to uncouple the displacement vector field and the rotation vector one (both of Lagrange type but displacement of second degree and rotation of first degree) and impose the Navier-Bernoulli condition by means of lagrange multipliers which in fact represent shear force field over the whole structure. This strategy worked. I put in attachment a script to perform stretching or bending of a 3D spring of circular section. I put also a mesh consistent with data (number of coils, height, radius...) entered in the py file. Ideally the mesh could easily be reconstructed by means of gmsh.py to be able to change the number of coils, height... but at that time I used pygmsh which does not work anymore on my computer. I also attach the vtk solutions for stretching and bending.

However, another point which was rather painful if I remember well was the difficulty to recover gmsh physical points. I can't remember if it is still the case in GetFEM since I haven't used it for some time now. In the spring example, I wrote a small function to find the dof node closest to a given point which works but is far being optimal. There may be a better way to do this...

Best regards
Jean-François Barthélémy


Le lun. 28 oct. 2024 à 08:59, Konstantinos Poulios <logari81@googlemail.com> a écrit :
Dear Jean-François and Yves,

The question about beam element modelling has appeared again in the getfem discourse forum:

In my understanding, there is no easy way to create a model of classical beam elements oriented in arbitrary directions in 2D/3D, with linear shape functions for the longitudinal displacement and Hermite shape functions for the transverse displacements. Even if assembling local stiffness matrices on discontinuous elements could be done, GetFEM has no simple mechanism of transforming local to global dofs and coupling the elements at nodes between elements with tangent discontinuities, without manually filling in rather complicated extension/reduction matrices for the mesh_fem.

The alternative is to use full Hermite shape functions for all deformation modes. I took over Jean-Francois'  partial implementations and adapted/extended it to the demo code posted here:
However I do not really get correct, mesh-independent, results. Maybe you have some quick hint on that.

Best regards
Kostas

On Tue, Nov 21, 2017 at 4:27 PM Jean-François Barthélémy <jfrancois.barthelemy@gmail.com> wrote:
Dear Yves,

Thank you very much. Now it works perfectly. It was indeed the problem of dimension of element_K which was at the origin of the error. I have downloaded the fixed version since the use of ":" did not work because of mismatch between dimensions (2x2 : 2x1).

The _expression_ "Grad_u:Grad_Test_u" is not equivalent concerning the contribution of the longitudinal part of the work since it also includes a part of work of the orthogonal displacement. The latter does not contribute to the extension work at first order but it rather contributes to the term related to the flexion expressed by means of its second derivative obtained by the Hessian. I've also built it and it works well (the analytical solution of the simple bending beam is retrieved).

I have been initially misled by the fact that the gradient of a vector for a 1D structure is actually a 2x2 in 2D (or 3x3 in 3D) matrix equal to du/ds@t where t is the vector "Normalized(element_K)" in GetFEM instead of just the vector du/ds as I would have expected since s is the only local coordinate. And the hessian is the third order tensor d2u/ds2@t@t instead of the vector d2u/ds2. But these vectors du/ds and d2u/ds2 can be obtained by contracting Grad_u with "Normalized(element_K)" on the one hand and Hess_u with "Normalized(element_K)@Normalized(element_K)" on the other hand. Once these operations are performed, everything is OK now and GetFEM works perfectly for arbitrary curvilinear structures.

Thank you again.

Best regards
Jean-François


2017-11-21 10:12 GMT+01:00 Yves Renard <Yves.Renard@insa-lyon.fr>:

Dear Jean-François,

The problem with the _expression_
"0.5*sqr(Normalized(element_K).(Grad_u*Normalized(element_K)))"

was that element_K is in fact declared as a matrix ( a 2x1 matrix here) and the dot product is restricted to vectors. I commited a fix to transform element_K as a vector for 1D elements. You can download this fix or use ":" instead of "." for the moment.

In fact, in that case, I think taht your _expression_ is equivalent to "Grad_u:Grad_Test_u" since the gradient is along the tangent.

Just a remark. It is of course allowed to use potential (it will be derived twice to obtain the tangent matrix), but I recommend to gives the weak formulation _expression_ not to have some surprise, especially for coupled problems (for instance if a coefficient is depending on another variable).

Best regards,

Yves



Le 21/11/2017 à 00:41, Jean-François Barthélémy a écrit :
Dear Yves,

Here are some small examples in Python on a simple beam [AB] made of 10 elements considered as a plane structure (problem in 2D). The mesh is loaded in GetFEM as a 2D mesh.

The file beam_bendmanual.py gives the solution of a beam in flexion and tension under the Navier-Bernoulli hypothesis. The formulation is built here by means of explicit indices in the assembly strings, which is what I would actually like to avoid with a complex structure made of arbitrary orientated beams. In this file, the displacement vector as well as the rotation of point A are blocked and a force loads the beam at point B. The stiffness are set to 1 for the sake of simplicity.

The file beam_bendgeneric.py gives the beggining of the attempt of a generic way to deal with the same problem (exploiting then the local direction). This file raises the exception mentioned in my previous email when trying to assembly. I haven't added the term of bending energy (involving the Hessian) nor the boundary condition blocking the rotation yet since they may not be as simple. For example the rotation should be blocked by imposing Grad_u.n=0 where n is the vector orthogonal to the beam direction but I am not sure that such a calculation of vectorial product 'n=ez x Normalized(element_K)' to get the vector n is available in the high-level assemblage syntax. Is it?

Thank you for your help.

Best regards
Jean-François


2017-11-20 15:00 GMT+01:00 Yves Renard <Yves.Renard@insa-lyon.fr>:

This should work, yes. If you can send me a small program that I can test, I can have a look to this problem.

Best regards,

Yves.



Le 20/11/2017 à 14:54, Jean-François Barthélémy a écrit :
Dear Yves,

Thank you. This is precisely the formulation I used but it raises the following problem (in python)

0.5*sqr(Normalized(element_K).(Grad_u*Normalized(element_K)))
------------------------------------------^
The second argument of the dot product has to be a vector.
logic_error exception caught
...
RuntimeError: (Getfem::InterfaceError) -- Error in getfem_generic_assembly.cc, line 8949 void getfem::ga_node_analysis(const string&, getfem::ga_tree&, const getfem::ga_workspace&, getfem::pga_tree_node, bgeot::size_type, bgeot::size_type, bool, bool, int): 
Error in assembly string

following a call such as
md.add_linear_generic_assembly_brick(mim,"0.5*sqr(Normalized(element_K).(Grad_u*Normalized(element_K)))")

Did I miss something?

I am sorry to bother you again.

Thanks

Best regards
Jean-François



2017-11-20 14:10 GMT+01:00 Yves Renard <Yves.Renard@insa-lyon.fr>:

Dear Jean-François,

For a vector variable 'u', each line of 'Grad_u' is the gradient of the ith component of 'u', each of them is tangent to the curve and length being the derivative with respect to the curvilinear abscissa. The linearized deformation is a priori  Normalized(element_K).(Grad_u * Normalized(element_K))


The formulas used to compute the gradient and the Hessian can be found here:

http://getfem.org/project/femdesc.html#geometric-transformations

http://getfem.org/project/appendixA.html#derivative-computation

The hessian of a vector valued variable is also the hessian of each component.

Best regards,

Yves.





Le 20/11/2017 à 02:10, Jean-François Barthélémy a écrit :
Dear Yves,

Thank you very much for your answer.

It's OK for scalar variables but I do not really understand how Grad_u is built when u is a displacement vector field of 3 components. I thought Grad_u would represent the vector du/ds ([dux/ds,duy/ds,duz/ds]) with s the local curvilinear abscissa (so that Normalized(element_K).Grad_u would give the linearized longitudinal deformation) but it seems that Grad_u is actually a 3x3 matrix field. Then I do not see how to build the longitudinal deformation. What would be the best way please? And by the way, what would be the right syntax to get the second derivative of the transverse displacement by means of Hermite elements and the Hessian?

Thank you again for your help.

Best regards
Jean-François



2017-11-17 20:52 GMT+01:00 Yves Renard <yves.renard@insa-lyon.fr>:

Dear Jean-François,

There is no specific tool yet for that.
You can have access to the tangent with 'element_K' in the generic assembly language (the unit tangent is then 'Normalized(element_K)')
If you define a scalar quantity "u" on your 1D structure, then "Grad_u" will be the gradient of the quantity in the sense that it is a tangent vector whose norm is the derivative of the qunatity along the curve. So that "Grad_u.Grad_Test_u" is still the stiffness term for a curvilinear second derivative. For a vector quantity "u", "Grad_u" is the componentwise gradient.

Best regard,

Yves.



----- Original Message -----
From: "Jean-François Barthélémy" <jfrancois.barthelemy@gmail.com>
To: getfem-users@nongnu.org
Sent: Friday, November 17, 2017 6:17:13 PM
Subject: [Getfem-users] Curvilinear structures in Getfem

Dear Getfem users,

I wonder whether it is possible to model simple linear elastic curvilinear
structures submitted to traction, bending, torsion etc... in 2D or 3D in
Getfem. I haven't found a way to have access to the tangential or normal
parts of vectors in the local basis of a beam and their derivatives with
respect to the curvilinear abscissa needed to build the formulation. Does
someone have an answer please?

Thanks in advance

Best regards
Jean-François


-- 

  Yves Renard (Yves.Renard@insa-lyon.fr)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------



-- 

  Yves Renard (Yves.Renard@insa-lyon.fr)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------



-- 

  Yves Renard (Yves.Renard@insa-lyon.fr)       tel : (33) 04.72.43.87.08
  Pole de Mathematiques, INSA-Lyon             fax : (33) 04.72.43.85.29
  20, rue Albert Einstein
  69621 Villeurbanne Cedex, FRANCE
  http://math.univ-lyon1.fr/~renard

---------

Attachment: spring.py
Description: Binary data

Attachment: spring.msh
Description: Binary data

Attachment: spring_bend.vtk
Description: Binary data

Attachment: spring_stretch.vtk
Description: Binary data

reply via email to
[Prev in Thread] Current Thread [Next in Thread]