 ## Solid and Structural Mechanics Demo 4: A planar beam

This code solves an elastic beam problem in structural mechanics theory. The libray contains the class Beam3DL2 that handles elastic beams in the space with 6 degrees of freedom for each node. These degrees of freedom are supported by the vertices of each element line.

• We start by defining an instance of class to deal with project data. We extract useful information from this instance: A verbosity and an output flag to deal with output levels. We next construct the Mesh and the Prescription instances.

 ``` IPF data("beam - 1.0",argv); int verbose = data.getVerbose(); int output_flag = data.getOutput(); Mesh ms(data.getMeshFile()); Prescription p(ms,data.getDataFile()); ```

• The linear system uses a symmetric matrix. For this, we we use an instance of class SkSMatrix<double>. Then we declare a vector that will contain the right-hand side and then the solution. Next, we read boundary conditions and prescribed concentrated loads from the prescription file. These are stored in vectors bc and bf respectively:

 ``` SkSMatrix A(ms); Vect b(ms), bc(ms), bf(ms); p.get(BOUNDARY_CONDITION,bc,0); p.get(POINT_FORCE,bf); ```

• The global matrix and right-hand side are, as usual, filled by looping over elements. For each element, an instance of class Beam3DL2 is considered. This class deals with quadrilateral elements. The member function Stiffness adds the contribution of the element stiffness matrix. The other member functions have already been considered in previous lessons. Note that the load vector is a copy of the concentrated load vector. No need for assembly here.

 ``` MeshElements(ms) { Beam3DL2 eq(theElement,0.1,0.1,0.1); eq.Stiffness(); eq.ElementAssembly(A); eq.ElementAssembly(b); } b = bf; ```

• We can now prescribe boundary conditions by penalty method and solve the resulting linear system by factorization and backsubstitution.

 ``` A.Prescribe(b,bc); A.solve(b); ```

• Let us now extract some additional information from the class. We want to output, for each element, the axial force, the bending moment, the twisting moment and the shear force. This is done simply by invoking corresponding member functions in the class Beam3DL2. Note that we need, before that, to construct the vector of node displacements. A dedicated constructor is made for that. Then we recall the member functions we need for post-processing.

 ``` Vect u(ms,3); Beam3DL2 eq(ms,uf,u); cout << "Element Axial Force Bending Moment Twisting Moment Shear Force" << endl; MeshElements(ms) { Beam3DL2 eq(theElement,0.1,0.1,0.1,b); cout << setw(8) << theElementLabel; cout << " " << eq.AxialForce(); cout << " " << eq.BendingMoment(); cout << " " << eq.TwistingMoment(); cout << " " << eq.ShearForce() << endl; } ```

• We finally show here how to create a deformed mesh for a graphical presentation. We have for this the function DeformMesh that modifies node coordinates by adding the displacement vector to the coordinates. The function can also multiply the displacement by an amplification factor which is equal here to 1 by default. Then we store the resulting mesh into a new file.

 ``` DeformMesh(ms,u); ms.put("deformed_beam.m"); ```

## An example

Let us run this program with an example of a straight bar [0,1]. We define the project file:

 ``` <date /> <author /> </info> <Project name="beam"> <mesh_file value="beam.dat"/> <data_file value="beam.dat" /> <verbose value="1" /> <output value="0"> <save value="1"/> <plot value="1"/> <init value="1"/> <plot_file value="beam.d"/> </Project> <Mesh dim="3" nb_dof="6"> <Nodes> 0.00000000e+00 0.00000000e+00 0.00000000e+00 111111 1.25000000e-01 0.00000000e+00 0.00000000e+00 0 2.50000000e-01 0.00000000e+00 0.00000000e+00 0 3.75000000e-01 0.00000000e+00 0.00000000e+00 0 5.00000000e-01 0.00000000e+00 0.00000000e+00 0 6.25000000e-01 0.00000000e+00 0.00000000e+00 0 7.50000000e-01 0.00000000e+00 0.00000000e+00 0 8.75000000e-01 0.00000000e+00 0.00000000e+00 0 1.00000000e+00 0.00000000e+00 0.00000000e+00 0 </Nodes> <Elements shape="line" nodes="2"> 1 2 1 2 3 1 3 4 1 4 5 1 5 6 1 6 7 1 7 8 1 8 9 1 </Elements> </Mesh> <Prescription> <PointForce x="1" dof="2">-0.01</PointForce> <PointForce x="1" dof="3"> 0.01</PointForce> </Prescription> </OFELI_File>```