Demo: Class TimeStepping: Solution of an Elastodynamics problem

We present a program to solve the time dependent linear elasticity equations (i.e. Elastodynamics) using the P1 finite element method for space discretization and the Newmark scheme for time discretization.

  • We start, as usual, by including the main header file of the library, and the header for Solid Mechanics problems. The program has, as argument the mesh file and the time step. We construct a Mesh instance with the flag that enables eliminating imposed boundary conditions. Note that the final time value is also set by using the OFELI global variable theFinalTime. Next, we declare the vectors to use in the calculations: u and v will contain the initial conditions for displacement and velocity, and u will contain at each time step the corresponding displacement vector. The initial displacement is set to the constant vector (0,-1)T, and the initial velocity to zero by default. 

    #include "OFELI.h"
#include "Solid.h"
using namespace OFELI;
                          
int main(int argc, char *argv[])
{
   theFinalTime = 1.;
   Mesh ms(argv[1],true);
   theTimeStep = atof(argv[2]);
                          
   Vect<double> u(ms), v(ms), bc(ms);
   u.set("-0.1",2);

  • Now we instantiate the equation class Elas2DT3 (Linearized planar elasticity using the P1 triangle). We also say that we the used terms in the equation are lumped mass, stiffness term (deviatoric and dialational). No body force is applied. We transmit all these data to the TimeStepping instance and choose as linear solver the conjugate gradient, preconditioned with the diagonal ILU.
    Since we intend to store for plotting displacement at all time steps, we store the initial displacements in a vtk file. 

       Elas2DT3 eq(ms);
   eq.setTerms(LUMPED_MASS|DEVIATORIC|DILATATION);
   TimeStepping ts(NEWMARK,theTimeStep,theFinalTime);
   ts.setPDE(eq);
   ts.setInitial(u,v);
   ts.setLinearSolver(CG_SOLVER,DILU_PREC);
   saveField(u,ms,"u-0.vtk",VTK);

  • We can now loop over time steps. For each one, we set boundary conditions (Here 0). We run the time step and store the resulted displacement vector in a vtk file. We can also store a deformed mesh in another file (More used in Solid Mechanics). 

       TimeLoop {
      bc.setTime(theTime);
      ts.setBC(bc);
      ts.runOneTimeStep();
      saveField(u,ms,"u-"+itos(theStep)+".vtk",VTK);
      Mesh dm(ms);
      DeformMesh(dm,u,1.);
      dm.save("mm-"+itos(theStep)+".m");
   }
   return 0;
}