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.getNbDOF());

• 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.getNbDOF()); p.get(BOUNDARY_CONDITION,bc,0); Vect<double> body_f(ms.getNbDOF()); p.get(BODY_FORCE,body_f); NodeVect<double> bound_f(ms,2); 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(ms,b,bc); a.Factor(); a.Solve(b);

• We can save the resulting displacement vector for plotting or any other purpose. For this, the class IOField is used.
  Note that we need to convert to class NodeVect<double> in order to save the vector with some mesh information.

  NodeVect<double> uf(ms,b,"Displacement"); if (save_flag) { IOField pl_file(data.getPlotFile(),OUT); pl_file.put(uf); }

• 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 ElementVect<double> instance will contain elementwise Von Mises stresses.

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

An example