Solid and Structural Mechanics Demo 2

Solid and Structural Mechanics Demo 2: A 3-D linear elasticity problem

This code solves three-dimensional linearized elastostatics problems. Its main structure is exactly the same as the 2-D code.

We start by defining an instance of class IPF to deal with project data. We extract useful information from this instance: An output flag to deal with output level and a save flag to say if the solution is to be saved in file. We then construct the Mesh and the Prescription instances. Note that the second argument in the mesh constructor means that imposed boundary conditions are to be removed from the list of unknowns. This is important in view of using iterative methods.

   IPF data("lelas3d - 1.2",argv[1]);
   int verbose = data.getVerbose(), output_flag = data.getOutput(), save_flag = data.getSave();
   Mesh ms(data.getMeshFile(),true);
   int nb_dof = ms.getDim();
   Prescription p(ms,data.getDataFile());

The linear system uses a symmetric matrix. For this, we we use an instance of class SpMatrix<double>. Then we declare a vector b and u that will contain respectively the right-hand side and the solution:
SpMatrix<double> A(ms); Vect<double> b(ms.getNbEq()), u(ms.getNbEq());
We then declare boundary conditions and loads. This is done by using class Prescription. Dirichlet boundary conditions are stored in vector bc, and body forces (ditributed loads) are in vector body_f. All these vectors can then by initialized from data as read in class Prescription.
Vect<double> bc(ms), body_f(ms), bound_f(ms); p.get(BOUNDARY_CONDITION,bc,0); p.get(BODY_FORCE,body_f,0); p.get(BOUNDARY_FORCE,bound_f,0);

The global matrix and right-hand side are filled by looping over elements. For each element, an instance of class Elas3DT4 is considered. This class deals with tetrahedral elements. Deviatoric and dilatational terms are constructed separately. The other member functions have already been considered in previous lessons. Note that we use the member function updateBC that takes into account boundary conditions at the element level.

   MeshElements(ms) {
      Elas3DT4 eq(theElement);
      eq.Deviator();
      eq.Dilatation();
      eq.BodyRHS(body_f);
      eq.updateBC(BCVect<double>(bc));
      eq.ElementAssembly(A);
      eq.ElementAssembly(b);
   }

We now solve the linear system by using the Conjugate Gradient method using the diagonal preconditioner.
CG(A,Prec<double>(A,DIAG_PREC),b,u,1000,toler,1);
and then construct the solution vector containing boundary conditions:
Vect<double> uf(ms); uf.insertBC(u,bc); if (output_flag > 0) cout << uf;
We save the resulting displacement vector for plotting or any other purpose. For this, the class IOField is used.
It is also classical in solid mechanics to represent the deformed domain. Here the function DeformMesh enables doing this. We then save the deformed mesh.
if (save_flag) { IOField pf(data.getPlotFile(),IOField::OUT); pf.put(uf); DeformMesh(ms,uf,1); ms.put(data.getProject()+"-1.m"); }

An example

Let us run this program with the data presented in the following project file:

<?xml version="1.0" encoding="ISO-8859-1" ?>
<OFELI_File>
<info>
   <title></title>
   <date></date>
   <author></author>
</info>
<Project name="beam">
   <verbose>1</verbose>
   <output>0</output>
   <save>1</save>
   <plot>1</plot>
   <init>1</init>
   <mesh_file>beam.m</mesh_file>
   <plot_file>beam.d</plot_file>
</Project>
<Prescription>
   <BodyForce dof="3">-1.0</BodyForce>
</Prescription>
</OFELI_File>

We note here that the file gives project parameters and prescriptions as well. Indeed, we have prescribed a vertical body force equal to -1 to the beam.