CxxFEM

The software which is named "fossils" here is based on a small-strain finite element code built with the help of the Gmsh SDK. This code is able to solve a wide range of linear statics problems. In practise, the code can be used through a python module named cxxfem.

Two examples of mechanical problems are given in the models/others folder of the installed version of fossils:

• beam2d.py: bending of a beam (2D plane stress),
• beam3d.py: combined bending and torsion of a beam (3D).

Let us have a look at beam2d.py.

The file starts with the importation of the cxxfem module as well as gmsh which is used for pre- and post-processing.

import cxxfem as fem
import gmsh

These two lines are folowed by a function named build_square which creates the problem geometry using a series Gmsh calls. The syntax is described in the Gmsh manual. It is also possible to load a .geo script from Gmsh (see beam3d.py which loads parallelepiped.geo).

The result of this function is a rectangle with several geometrical entities named with explicit names ("top" is the top curve, "interior" is the interior surface of the beam, etc.). These names are also called "physical groups" in Gmsh. They are used to assign properties and boundary conditions to a list of geometrical entities.

The function is followed by the __main__ part of the script, which is called when the script is executed by python.

Every problem definition starts with the creation of a Problem object.

pbl = fem.Problem()

The object is stored in a python variable named pbl. The fem. prefix refers to the cxxfem module which was imported as fem at the beginning of the file. This object is the main object containing the main data related to the finite element analysis.

Then, the build_square function described ealier is called with particular parameters (the beam size and the mesh size)

build_square(0.0, 0.0, 10.0, 1.0, 80, 8, True, True)

The next step is to define material properties. This is done by adding a Medium object to the Problem. The target physical group ("interior") is also given, as well as the values of the Young's modulus and the Poisson coefficient of the linear elastic model:

m = fem.Medium(pbl, "interior", 100.0, 0.3)

Next, we can define boundary conditions:

d1 = fem.Dirichlet(pbl, "left", "X", 0.0)
d2 = fem.Dirichlet(pbl, "left", "Y", 0.0)
n3 = fem.Neumann(pbl, "top", "Y", -1e-3)
d4 = fem.Dirichlet(pbl, "interior", "Z", 0.0)

Their type can be either Dirichlet (the value of the unknown is prescribed) or Neumann (the derivative of the unknown is prescribed).

The syntax is rather simple. For example, the first line means that we want to prescribe the value of the X component of the displacement to be 0 and the "left" side of the beam. Similarly, the Y component of the displacement is fixed to 0 too, so that the left side is fully clamped. The Neumann condition means that we want to apply a uniformly distributed load of -1e-3 along Y on the top side of the beam. Finally, the Z component of the dispacement field is fixed to 0 since it is a 2D problem.

At this stage, the Problem is fully defined. It should be solved by a Solver object.

The solver is defined by the following line. It takes the target problem as argument:

solver = fem.Solver(pbl)

The solver is then executed:

solver.solve()

This command starts the calculations. When they are finished, the results can be displayed and analysed. This part is managed by a Post object, which takes the target solver as argument:

post = fem.Post(solver)

The post-processing views (stress/strain fields calculation from the displacements previously computed) are built with the following command:

post.build_views()

The views that will be displayed later in Gmsh can be selected with:

post.show_views( [ "stress_tensor", "force_vector" ] )

It is also possible (but not mandatory) to write the results to disk in the current workspace folder:

post.write()

Since the FE code solves small-strain problems, the deformation of the mesh should be magnified in order to be visible. This is done with the following command:

post.deform(5)  

The Post object is also able to extract results on given physical groups. For example, we can extract the values of the force_vector on the nodes lying of the left side of the beam:

forces = post.probe("force_vector", "left")

The values are returned in a python list such as [v1x, v1y, v1z, v2x, v2y, v2z, ..., vnx, vny, vnZ] where [v1x, v1y, v1z] are the values of the 3 components of the force vector at the first node of the physical group; [v1x, v1y, v1z] are the values related to the second node, and so on...

It is thus easy to calculate the resultant vector whith the following python commands:

print(f'Rx = {sum(forces[::3])} N')
print(f'Ry = {sum(forces[1::3])} N')

Similar commands can be used to extract the vertical displacement of the beam:

disp = post.probe("Y", "pt2")
print(f'tip displacement = {disp[0]} mm')

The last command of the script starts Gmsh and displays the results:

post.view()