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 <[email protected]>: > > 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 <[email protected]>: > >> >> 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 <[email protected]>: >> >>> >>> 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" <[email protected]> >>> To: [email protected] >>> 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 ([email protected]) 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 ([email protected]) 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 > > --------- > >
#!/usr/bin/env python
from numpy import *
from getfem import *
DIRICHLET=1;NEUMANN=2
m=Mesh('import', 'gmsh', 'beam.msh')
mfu = MeshFem(m,2);mfu.set_fem(Fem('FEM_HERMITE(1)'))
mfl = MeshFem(m,1);mfl.set_fem(Fem('FEM_PK(1,1)'))
mim = MeshIm(m,2)
md=Model('real')
md.add_fem_variable('u',mfu)
md.add_initialized_data('F',[1.,1.])
md.add_linear_generic_assembly_brick(mim,"0.5*sqr(Grad_u(1,1))")
md.add_linear_generic_assembly_brick(mim,"0.5*sqr(Hess_u(2,1,1))")
md.add_source_term_brick(mim,'u','F',NEUMANN)
md.add_Dirichlet_condition_with_multipliers(mim,'u',mfu,DIRICHLET)
md.add_filtered_fem_variable('lag',mfl,DIRICHLET)
md.add_linear_generic_assembly_brick(mim,"lag*Grad_u(2,1)",DIRICHLET)
md.solve()
mfu.export_to_vtk('beam.vtk','ascii',mfu,md.variable('u'),'u')
beam.msh
Description: Mesh model
#!/usr/bin/env python
from numpy import *
from getfem import *
DIRICHLET=1;NEUMANN=2
m=Mesh('import', 'gmsh', 'beam.msh')
m.display()
mfu = MeshFem(m,2);mfu.set_fem(Fem('FEM_HERMITE(1)'))
mfl = MeshFem(m,1);mfl.set_fem(Fem('FEM_PK(1,1)'))
mim = MeshIm(m,2)
md=Model('real')
md.add_fem_variable('u',mfu)
md.add_initialized_data('F',[1.,1.])
md.add_linear_generic_assembly_brick(mim,"0.5*sqr(Normalized(element_K).(Grad_u*Normalized(element_K)))")
# + bending energy (with Hessian) = int(0.5*sqr(d2(u_ortho)/ds2))
md.add_source_term_brick(mim,'u','F',NEUMANN)
md.add_Dirichlet_condition_with_multipliers(mim,'u',mfu,DIRICHLET)
# + boundary condition to block the rotation ie d(u_ortho)/ds
md.solve()
# mfu.export_to_vtk('beam.vtk','ascii',mfu,md.variable('u'),'u')
