Solid and Structural Mechanics Demo 1

Solid and Structural Mechanics Demo 1: A linear elasticity problem with planar deformations

This code solves linearized elastostatics problems with planar deformations.

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.
IPF data("lelas2d - 1.1 ",argv[1]); int output_flag = data.getOutput(); int save_flag = data.getSave(); 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:
SkSMatrix<double> A(ms); Vect<double> b(ms);
We can then declare boundary conditions and loads. This is done by using class Prescription. Dirichlet boundary conditions are stored in vector bc, body forces (ditributed loads) are in vector body_f and boundary forces (tractions or pressures on boundary) are in vector bound_f. Note that we have chosen to store boundary forces nodewise. All these vectors can then by initialized from data as read in class Prescription.
Vect<double> bc(ms); p.get(BOUNDARY_CONDITION,bc,0); Vect<double> body_f(ms); p.get(BODY_FORCE,body_f); Vect<double> bound_f(ms); p.get(TRACTION,bound_f,0);

The global matrix and right-hand side are, as usual, filled by looping over elements. For each element, an instance of class Elas2DQ4 is considered. This class deals with quadrilateral elements. Deviatoric and dilatational terms are constructed separately. The other member functions have already been considered in previous lessons.

   MeshElements(ms) {
      Elas2DQ4 eq(theElement);
      eq.Deviator();
      eq.Dilatation();
      eq.BodyRHS(Vect<double>(theElement,body_f));
      eq.ElementAssembly(A);
      eq.ElementAssembly(b);
   }

Since we want to take account for boundary forces, we have to loop over sides (here edges) and construct element loads. This has also been considered in previous lessons.
MeshSides(ms) { Elas2DQ4 eq(theSide); eq.BoundaryRHS(Vect<double>(theSide,bound_f)); eq.SideAssembly(b); }
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);
We can save the resulting displacement vector for plotting or any other purpose. For this, the class IOField is used.
if (save_flag) { IOField pl_file(data.getPlotFile(),IOField::OUT); pl_file.put(b); }

Let us illustrate how to do some postprocessing. We want here to compute principal and Von Mises stresses. For this, we have to use the member function Stress. The element principal stresses are stored in the instance ste of class LocalVect<double,3> while an Vect<double> instance will contain elementwise Von Mises stresses.

   Vect<double> st(ms,"Principal Stress",0.0,3,ELEMENT_FIELD);
   Vect<double> vm(ms,"Von-Mises Stress",0.0,3,ELEMENT_FIELD);
   MeshElements(ms) {
      LocalVect<double,3> ste;
      Elas2DQ4 eq(theElement,b);
      eq.Stress(ste,vm(theElement->getLabel()));
      st(theElementLabel,1) = ste(1);
      st(theElementLabel,2) = ste(2);
      st(theElementLabel,3) = ste(3);
   }

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>1</output>
   <save>1</save>
   <plot>1</plot>
   <init>1</init>
   <mesh_file>beam.m</mesh_file>
   <plot_file>beam.d</plot_file>
   <prescription_file>beam.pr</prescription_file>
</Project>
</OFELI_File>

The only thing that may differ from other seen project files is that prescription data is given in a separate file: beam.pr. Let us give this one:

<?xml version="1.0" encoding="ISO-8859-1" ?>
<OFELI_File>
<info>
   <title></title>
   <date></date>
   <author></author>
</info>
<Prescription>
   <BoundaryCondition code="1">0.0</BoundaryCondition>
   <Traction code="1" dof="2">-1000.0</Traction>
</Prescription>
</OFELI_File>

Clearly, we have prescribed by this file a (Dirichlet) boundary condition that enforces a value of 0.0 for degrees of freedom with code 1. In addition, we have prescribed a traction of value -1000 for the second degrees of freedom (y-traction) having code 1.