Beam Model Simulation in GETFEM

Questions and issues
6

it is possible to define beam elements using Hermite shape functions for the bending part and linear element for the axial deformations. These two components will give you the same local element stiffness matrix as in the usual beam theory. If the beam elements are in arbitrary orientations in 2D or 3D space you will need to do the local to global coordinate transformation and the coupling between the elements manually.

There is a previous discussion about beam elements here
https://lists.gnu.org/archive/html/getfem-users/2017-11/msg00026.html

Based on Jean Francois’ code from this reference, defining beam elements using only Hermite shape functions would look something like this (but still the code does not yield correct results, needs some debugging)

#!/usr/bin/env python

import getfem as gf
RG_FIXED = 99;
m = gf.Mesh("import", "structured",
               f"GT='GT_PK(1,1)';ORG=[0.];SIZES=[5.];NSUBDIV=[10]")
#m.transform([[1.],[0.5]]) # inclined beam
m.transform([[1.],[0.]]) # horizontal beam
m.set_region(RG_FIXED, m.outer_faces_in_box([-1e-5,-1e-5], [1e-5,1e-5], 1)) # select lower left corner
m.export_to_vtu("mesh.vtu")
mfu = gf.MeshFem(m,2)
mfu.set_fem(gf.Fem("FEM_HERMITE(1)"))
mfl = gf.MeshFem(m,1)
mfl.set_fem(gf.Fem("FEM_PK(1,1)"))
mim = gf.MeshIm(m,2)
mimd = gf.MeshImData(mim,-1,2) # "-1" means the entire domain
md = gf.Model("real")
md.add_fem_variable("u",mfu)
md.add_im_data("Tvec", mimd)
md.add_im_data("Nvec", mimd)
md.add_macro("Tvecdef", "Normalized(element_K)")
md.add_macro("Nvecdef", "[0,-1;1,0]*Normalized(element_K)")
md.set_variable("Tvec", md.interpolation("Tvecdef", mimd))
md.set_variable("Nvec", md.interpolation("Nvecdef", mimd))
md.add_linear_term(mim,"0.5*sqr((Grad(u)*Tvec)).Tvec")
md.add_linear_term(mim,"0.5*sqr(((Hess(u).Tvec).Tvec).Nvec)")
md.add_linear_term(mim,"0.01*Test_u(2)") # uniform vertical load
md.add_Dirichlet_condition_with_multipliers(mim,"u",mfu,RG_FIXED)
#md.add_linear_term(mim,"1e6*u.Test_u", RG_FIXED) # fixed position left node by penalization
md.add_filtered_fem_variable("lag_slope", mfl, RG_FIXED)
md.add_linear_term(mim,"lag_slope*Nvecdef.Grad(u).Tvecdef", RG_FIXED)
md.solve()
mfu.export_to_vtu("beam.vtu","ascii",mfu,md.variable("u"),"u")