Example codes
and tutorials

The (Not-so) Quick Guide

List of tutorials/demo codes

Single-Physics Problems

Poisson

Adaptivity illustrated for Poisson

Advection-Diffusion

Unsteady heat equation

Linear wave equation

The Young-Laplace equation

Navier-Stokes

Free-surface Navier-Stokes

Axisymmetric Navier-Stokes

Solid mechanics

Beam structures

Shell structures

Multi-physics Problems

Fluid-structure interaction

Boussinesq convection

Steady thermoelasticity

Methods-based example codes and tutorials

Mesh generation

Linear solvers and preconditioners

Visualisation of the results

Parallel processing

How to write a new element

How to write a new refineable element

Default nonlinear solvers -- the sequence of action functions

...

Documentation

FE theory and top-down discussion of the data structure

The (Not-so) Quick Guide

Comprehensive bottom-up discussion of the data structure

List of available structured and unstructured meshes

Linear solvers and preconditioners

Visualisation of the results

Parallel processing

Coding conventions and C++ style

Creating documentation

Optimisation - robustness vs. "raw speed"

Linear vs. nonlinear problems

Storing shape functions

Changing the default "full" integration scheme

Disabling the ALE formulation of unsteady equations

C vs. C++ output

Different sparse assembly techniques and the STL memory pool

Publications

Publications

Talks

Journal publications

Theses

Picture show

Download

Copyright

Download/installation instructions

Download page

FAQ & Contact

FAQ

Change log

Bugs and other known problems

Completeness of the library & our "To-Do List"

Contact the developers

Get involved

 

---

Beta release! Please note that the library has not been "officially" released. While we continue to work on the documentation, these web pages are likely to contain broken links and documents in draft form. Please send an email to oomph-lib AT maths DOT man DOT ac DOT uk if you wish to be informed of the library's "official" release.

---

shock_disk.cc

Go to the documentation of this file.

//LIC// ====================================================================
//LIC// This file forms part of oomph-lib, the object-oriented, 
//LIC// multi-physics finite-element library, available 
//LIC// at http://www.oomph-lib.org.
//LIC// 
//LIC//           Version 0.90. August 3, 2009.
//LIC// 
//LIC// Copyright (C) 2006-2009 Matthias Heil and Andrew Hazel
//LIC// 
//LIC// This library is free software; you can redistribute it and/or
//LIC// modify it under the terms of the GNU Lesser General Public
//LIC// License as published by the Free Software Foundation; either
//LIC// version 2.1 of the License, or (at your option) any later version.
//LIC// 
//LIC// This library is distributed in the hope that it will be useful,
//LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
//LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
//LIC// Lesser General Public License for more details.
//LIC// 
//LIC// You should have received a copy of the GNU Lesser General Public
//LIC// License along with this library; if not, write to the Free Software
//LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
//LIC// 02110-1301  USA.
//LIC// 
//LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
//LIC// 
//LIC//====================================================================
// Driver for large-displacement elasto-dynamic test problem:
// Circular disk impulsively loaded by compressive load.

#include <iostream>
#include <fstream>
#include <cmath>

//My own includes
#include "generic.h"
#include "solid.h"

//Need to instantiate templated mesh
#include "meshes/quarter_circle_sector_mesh.h"

using namespace std;

using namespace oomph;
 
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////


//================================================================
/// Global variables
//================================================================
namespace Global_Physical_Variables
{
 /// Pointer to constitutive law
 ConstitutiveLaw* Constitutive_law_pt;

 /// Elastic modulus
 double E=1.0;

 /// Poisson's ratio
 double Nu=0.3;

 /// Uniform pressure
 double P = 0.00;

 /// Constant pressure load
 void constant_pressure(const Vector<double> &xi,const Vector<double> &x,
                        const Vector<double> &n, Vector<double> &traction)
 {
  unsigned dim = traction.size();
  for(unsigned i=0;i<dim;i++)
   {
    traction[i] = -P*n[i];
   }
 }

}



///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////



//================================================================
/// Elastic quarter circle sector mesh with functionality to
/// attach traction elements to the curved surface. We "upgrade"
/// the RefineableQuarterCircleSectorMesh to become an
/// SolidMesh and equate the Eulerian and Lagrangian coordinates,
/// thus making the domain represented by the mesh the stress-free 
/// configuration. 
/// \n\n
/// The member function \c make_traction_element_mesh() creates
/// a separate mesh of SolidTractionElements that are attached to the
/// mesh's curved boundary (boundary 1). 
//================================================================
template <class ELEMENT>
class ElasticRefineableQuarterCircleSectorMesh :
 public virtual RefineableQuarterCircleSectorMesh<ELEMENT>,
 public virtual SolidMesh
{


public:

 /// \short Constructor: Build mesh and copy Eulerian coords to Lagrangian
 /// ones so that the initial configuration is the stress-free one.
 ElasticRefineableQuarterCircleSectorMesh<ELEMENT>(GeomObject* wall_pt,
                                         const double& xi_lo,
                                         const double& fract_mid,
                                         const double& xi_hi,
                                         TimeStepper* time_stepper_pt=
                                         &Mesh::Default_TimeStepper) :
  RefineableQuarterCircleSectorMesh<ELEMENT>(wall_pt,xi_lo,fract_mid,xi_hi,
                                             time_stepper_pt)
  {
#ifdef PARANOID
   /// Check that the element type is derived from the SolidFiniteElement
   SolidFiniteElement* el_pt=dynamic_cast<SolidFiniteElement*>
    (finite_element_pt(0));
   if (el_pt==0)
    {
     throw OomphLibError(
      "Element needs to be derived from SolidFiniteElement\n",
      "ElasticRefineableQuarterCircleSectorMesh::constructor()",
      OOMPH_EXCEPTION_LOCATION);
    }
#endif

   // Make the current configuration the undeformed one by
   // setting the nodal Lagrangian coordinates to their current
   // Eulerian ones
   set_lagrangian_nodal_coordinates();
  }


 /// Function to create mesh made of traction elements
 void make_traction_element_mesh(SolidMesh*& traction_mesh_pt)
  {

   // Make new mesh
   traction_mesh_pt=new SolidMesh;

   // Loop over all elements on boundary 1:
   unsigned b=1;
   unsigned n_element = this->nboundary_element(b);
   for (unsigned e=0;e<n_element;e++)
    {
     // The element itself:
     FiniteElement* fe_pt = this->boundary_element_pt(b,e);
     
     // Find the index of the face of element e along boundary b
     int face_index = this->face_index_at_boundary(b,e);
     
     // Create new element
     traction_mesh_pt->add_element_pt(new SolidTractionElement<ELEMENT>
                                      (fe_pt,face_index));
    }
  }


 /// Function to wipe and re-create mesh made of traction elements
 void remake_traction_element_mesh(SolidMesh*& traction_mesh_pt)
  {

   // Wipe existing mesh (but don't call it's destructor as this
   // would wipe all the nodes too!)
   traction_mesh_pt->flush_element_and_node_storage();

   // Loop over all elements on boundary 1:
   unsigned b=1;
   unsigned n_element = this->nboundary_element(b);
   for (unsigned e=0;e<n_element;e++)
    {
     // The element itself:
     FiniteElement* fe_pt = this->boundary_element_pt(b,e);
     
     // Find the index of the face of element e along boundary b
     int face_index = this->face_index_at_boundary(b,e);
     
     // Create new element
     traction_mesh_pt->add_element_pt(new SolidTractionElement<ELEMENT>
                                      (fe_pt,face_index));
    }
  }


};





/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////



//====================================================================== 
/// "Shock" wave propagates through an impulsively loaded
/// circular disk.
//====================================================================== 
template<class ELEMENT, class TIMESTEPPER>
class DiskShockWaveProblem : public Problem
{

public:

 /// Constructor:
 DiskShockWaveProblem();

 /// \short Run the problem; specify case_number to label output
 /// directory
 void run(const unsigned& case_number);
 
 /// Access function for the solid mesh
 ElasticRefineableQuarterCircleSectorMesh<ELEMENT>*& solid_mesh_pt() 
  {
   return Solid_mesh_pt;  
  } 

 /// Access function for the mesh of surface traction elements
 SolidMesh*& traction_mesh_pt() 
  {
   return Traction_mesh_pt;  
  } 

 /// Doc the solution
 void doc_solution();

 /// Update function (empty)
 void actions_after_newton_solve() {}

 /// Update function (empty)
 void actions_before_newton_solve() {}

 /// \short Actions after adaption: Kill and then re-build the traction 
 /// elements on boundary 1 and re-assign the equation numbers
 void actions_after_adapt();

 /// \short Doc displacement and velocity: label file with before and after
 void doc_displ_and_veloc(const int& stage=0);

 /// \short Dump problem-specific parameters values, then dump
 /// generic problem data.
 void dump_it(ofstream& dump_file);

 /// \short Read problem-specific parameter values, then recover
 /// generic problem data.
 void restart(ifstream& restart_file);

private:

 // Output
 DocInfo Doc_info;

 /// Trace file
 ofstream Trace_file;
 
 /// Vector of pointers to nodes whose position we're tracing
 Vector<Node*> Trace_node_pt;

 /// Pointer to solid mesh
 ElasticRefineableQuarterCircleSectorMesh<ELEMENT>* Solid_mesh_pt;

 /// Pointer to mesh of traction elements
 SolidMesh* Traction_mesh_pt;

};





//====================================================================== 
/// Constructor
//====================================================================== 
template<class ELEMENT, class TIMESTEPPER>
DiskShockWaveProblem<ELEMENT,TIMESTEPPER>::DiskShockWaveProblem() 
{


 //Allocate the timestepper
 add_time_stepper_pt(new TIMESTEPPER);
 
 // Set coordinates and radius for the circle that defines 
 // the outer curvilinear boundary of the domain
 double x_c=0.0;
 double y_c=0.0;
 double r=1.0;
 
 // Build geometric object that specifies the fish back in the
 // undeformed configuration (basically a deep copy of the previous one)
 GeomObject* curved_boundary_pt=new Circle(x_c,y_c,r,time_stepper_pt());

 // The curved boundary of the mesh is defined by the geometric object
 // What follows are the start and end coordinates on the geometric object:
 double xi_lo=0.0;
 double xi_hi=2.0*atan(1.0);

 // Fraction along geometric object at which the radial dividing line
 // is placed
 double fract_mid=0.5;

 //Now create the mesh
 solid_mesh_pt() = new ElasticRefineableQuarterCircleSectorMesh<ELEMENT>(
  curved_boundary_pt,xi_lo,fract_mid,xi_hi,time_stepper_pt());

 // Set up trace nodes as the nodes on boundary 1 (=curved boundary) in 
 // the original mesh (they exist under any refinement!) 
 unsigned nnod0=solid_mesh_pt()->nboundary_node(0);
 unsigned nnod1=solid_mesh_pt()->nboundary_node(1);
 Trace_node_pt.resize(nnod0+nnod1);
 for (unsigned j=0;j<nnod0;j++)
  {
   Trace_node_pt[j]=solid_mesh_pt()->boundary_node_pt(0,j);
  }
 for (unsigned j=0;j<nnod1;j++)
  {
   Trace_node_pt[j+nnod0]=solid_mesh_pt()->boundary_node_pt(1,j);
  }

 // Build traction element mesh
 solid_mesh_pt()->make_traction_element_mesh(traction_mesh_pt());
 
 // Solid mesh is first sub-mesh
 add_sub_mesh(solid_mesh_pt());

 // Traction mesh is first sub-mesh
 add_sub_mesh(traction_mesh_pt());

 // Build combined "global" mesh
 build_global_mesh();


 // Create/set error estimator
 solid_mesh_pt()->spatial_error_estimator_pt()=new Z2ErrorEstimator;
  
 // Fiddle with adaptivity targets and doc
 solid_mesh_pt()->max_permitted_error()=0.006; //0.03;
 solid_mesh_pt()->min_permitted_error()=0.0006;// 0.0006; //0.003;
 solid_mesh_pt()->doc_adaptivity_targets(cout);

 // Pin the bottom in the vertical direction
 unsigned n_bottom = solid_mesh_pt()->nboundary_node(0);

 //Loop over the node
 for(unsigned i=0;i<n_bottom;i++)
  {
   solid_mesh_pt()->boundary_node_pt(0,i)->pin_position(1);
  }

 // Pin the left edge in the horizontal direction
 unsigned n_side = solid_mesh_pt()->nboundary_node(2);
 //Loop over the node
 for(unsigned i=0;i<n_side;i++)
  {
   solid_mesh_pt()->boundary_node_pt(2,i)->pin_position(0);
  }

 //Find number of elements in solid mesh
 unsigned n_element =solid_mesh_pt()->nelement();
  
 //Set parameters and function pointers for solid elements
 for(unsigned i=0;i<n_element;i++)
  {
   //Cast to a solid element
   ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));

   //Set the constitutive law
   el_pt->constitutive_law_pt() =
    Global_Physical_Variables::Constitutive_law_pt;
   
   // Switch on inertia
   el_pt->unsteady()=true;

   // Use MacroElement representation for 
   // Lagrangian coordinates of newly created 
   // nodes hierher simply retained to let self tests pass during
   // re-development of code
   el_pt->use_undeformed_macro_element_for_new_lagrangian_coords()
    =true;
  }

 // Pin the redundant solid pressures
 PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
  solid_mesh_pt()->element_pt());

 //Find number of elements in traction mesh
 n_element=traction_mesh_pt()->nelement();
  
 //Set function pointers for traction elements
 for(unsigned i=0;i<n_element;i++)
  {
   //Cast to a solid traction element
   SolidTractionElement<ELEMENT> *el_pt = 
    dynamic_cast<SolidTractionElement<ELEMENT>*>
    (traction_mesh_pt()->element_pt(i));

   //Set the traction function
   el_pt->traction_fct_pt() = Global_Physical_Variables::constant_pressure;
  }

 //Attach the boundary conditions to the mesh
 cout << assign_eqn_numbers() << std::endl; 

 // Refine uniformly
 refine_uniformly();
 refine_uniformly();
 refine_uniformly();
 

 // Now the non-pinned positions of the SolidNodes will have been
 // determined by interpolation. This is appropriate for uniform
 // refinements once the code is up and running since we can't place
 // new SolidNodes at the positions determined by the MacroElement.
 // However, here we want to update the nodes to fit the exact
 // intitial configuration.

 // Update all solid nodes based on the Mesh's Domain/MacroElement
 // representation
 bool update_all_solid_nodes=true;
 solid_mesh_pt()->node_update(update_all_solid_nodes);

 // Now set the Eulerian equal to the Lagrangian coordinates
 solid_mesh_pt()->set_lagrangian_nodal_coordinates();
 
}






//==================================================================
/// Kill and then re-build the traction elements on boundary 1,
/// pin redundant pressure dofs and re-assign the equation numbers.
//==================================================================
template<class ELEMENT, class TIMESTEPPER>
void DiskShockWaveProblem<ELEMENT,TIMESTEPPER>::actions_after_adapt()
{ 
 // Wipe and re-build traction element mesh
 solid_mesh_pt()->remake_traction_element_mesh(traction_mesh_pt());
 
 // Re-build combined "global" mesh
 rebuild_global_mesh();

 //Find number of elements in traction mesh
 unsigned n_element=traction_mesh_pt()->nelement();
  
 //Loop over the elements in the traction element mesh
 for(unsigned i=0;i<n_element;i++)
  {
   //Cast to a solid traction element
   SolidTractionElement<ELEMENT> *el_pt = 
    dynamic_cast<SolidTractionElement<ELEMENT>*>
    (traction_mesh_pt()->element_pt(i));

   //Set the traction function
   el_pt->traction_fct_pt() = Global_Physical_Variables::constant_pressure;
  }

 // Pin the redundant solid pressures
 PVDEquationsBase<2>::pin_redundant_nodal_solid_pressures(
  solid_mesh_pt()->element_pt());


 //Do equation numbering
 cout << assign_eqn_numbers() << std::endl; 

}


//==================================================================
/// Doc the solution
//==================================================================
template<class ELEMENT, class TIMESTEPPER>
void DiskShockWaveProblem<ELEMENT,TIMESTEPPER>::doc_solution()
{
 // Number of plot points
 unsigned npts;
 npts=5; 

 // Output shape of deformed body
 ofstream some_file;
 char filename[100];
 sprintf(filename,"%s/soln%i.dat",Doc_info.directory().c_str(),
         Doc_info.number());
 some_file.open(filename);
 solid_mesh_pt()->output(some_file,npts);
 some_file.close();


 // Output traction
 unsigned nel=traction_mesh_pt()->nelement();
 sprintf(filename,"%s/traction%i.dat",Doc_info.directory().c_str(),
         Doc_info.number());
 some_file.open(filename);
 Vector<double> unit_normal(2);
 Vector<double> traction(2);
 Vector<double> x_dummy(2);
 Vector<double> s_dummy(1);
 for (unsigned e=0;e<nel;e++)
  {
   some_file << "ZONE " << std::endl;
   for (unsigned i=0;i<npts;i++)
    {
     s_dummy[0]=-1.0+2.0*double(i)/double(npts-1);
     SolidTractionElement<ELEMENT>* el_pt=
      dynamic_cast<SolidTractionElement<ELEMENT>*>(
       traction_mesh_pt()->finite_element_pt(e));
     el_pt->outer_unit_normal(s_dummy,unit_normal);
     el_pt->traction(s_dummy,traction);
     el_pt->interpolated_x(s_dummy,x_dummy);
     some_file << x_dummy[0] << " " << x_dummy[1] << " " 
               << traction[0] << " " << traction[1] << " "  
               << std::endl;
    }
  }
 some_file.close(); 

 // Doc displacement and velocity
 doc_displ_and_veloc();

 // Get displacement as a function of the radial coordinate
 // along boundary 0
 {

  // Number of elements along boundary 0:
  unsigned nelem=solid_mesh_pt()->nboundary_element(0);

  // Open files
  sprintf(filename,"%s/displ%i.dat",Doc_info.directory().c_str(),
          Doc_info.number());
  some_file.open(filename);
  
  Vector<double> s(2);
  Vector<double> x(2);
  Vector<double> dxdt(2);
  Vector<double> xi(2);
  Vector<double> r_exact(2);
  Vector<double> v_exact(2);

  for (unsigned e=0;e<nelem;e++)
   {
    some_file << "ZONE " << std::endl;
    for (unsigned i=0;i<npts;i++)
     {
      // Move along bottom edge of element
      s[0]=-1.0+2.0*double(i)/double(npts-1);
      s[1]=-1.0;

      // Get pointer to element
      SolidFiniteElement* el_pt=dynamic_cast<SolidFiniteElement*>
       (solid_mesh_pt()->boundary_element_pt(0,e));
      
      // Get Lagrangian coordinate
      el_pt->interpolated_xi(s,xi);

      // Get Eulerian coordinate
      el_pt->interpolated_x(s,x);

      // Get velocity 
      el_pt->interpolated_dxdt(s,1,dxdt);
  
      // Plot radial distance and displacement
      some_file << xi[0] << " " << x[0]-xi[0] << " " 
                << dxdt[0] << std::endl;
     }
   }
  some_file.close(); 
 }

  
 // Write trace file
 Trace_file << time_pt()->time()  << " ";
 unsigned ntrace_node=Trace_node_pt.size();
 for (unsigned j=0;j<ntrace_node;j++)
  {
   Trace_file << sqrt(pow(Trace_node_pt[j]->x(0),2)+
                      pow(Trace_node_pt[j]->x(1),2)) << " ";
  }
 Trace_file << std::endl;

 // Output principal stress vectors at the centre of all elements
 SolidHelpers::doc_2D_principal_stress<ELEMENT>(Doc_info,solid_mesh_pt());

//  // Write restart file
//  sprintf(filename,"%s/restart%i.dat",Doc_info.directory().c_str(),
//          Doc_info.number());
//  some_file.open(filename);
//  dump_it(some_file);
//  some_file.close();
 

 cout << "Doced solution for step " 
      << Doc_info.number() 
      << std::endl << std::endl << std::endl;
}






//==================================================================
/// Doc displacement and veloc in displ_and_veloc*.dat.
/// The int stage defaults to 0, in which case the '*' in the
/// filename is simply the step number specified by the Problem's
/// DocInfo object. If it's +/-1, the word "before" and "after"
/// get inserted. This allows checking of the veloc/displacment
/// interpolation during adaptive mesh refinement.
//==================================================================
template<class ELEMENT, class TIMESTEPPER>
void DiskShockWaveProblem<ELEMENT,TIMESTEPPER>::doc_displ_and_veloc(
 const int& stage)
{

 ofstream some_file;
 char filename[100];

 // Number of plot points
 unsigned npts;
 npts=5; 

 // Open file
 if (stage==-1)
  {
   sprintf(filename,"%s/displ_and_veloc_before%i.dat",
           Doc_info.directory().c_str(),Doc_info.number());
  }
 else if (stage==1)
  {
   sprintf(filename,"%s/displ_and_veloc_after%i.dat",
           Doc_info.directory().c_str(),Doc_info.number());
  }
 else 
  {
   sprintf(filename,"%s/displ_and_veloc%i.dat",
           Doc_info.directory().c_str(),Doc_info.number());
  }
 some_file.open(filename);

 Vector<double> s(2),x(2),dxdt(2),xi(2),displ(2);

 //Loop over solid elements
 unsigned nel=solid_mesh_pt()->nelement();
 for (unsigned e=0;e<nel;e++)
  {
   some_file << "ZONE I=" << npts << ", J=" << npts << std::endl;
   for (unsigned i=0;i<npts;i++)
    {
     s[0]=-1.0+2.0*double(i)/double(npts-1);
     for (unsigned j=0;j<npts;j++)
      {
       s[1]=-1.0+2.0*double(j)/double(npts-1);

       // Cast to full element type
       ELEMENT* el_pt=dynamic_cast<ELEMENT*>(solid_mesh_pt()->
                                             finite_element_pt(e));

       // Eulerian coordinate
       el_pt->interpolated_x(s,x);

       // Lagrangian coordinate
       el_pt->interpolated_xi(s,xi);

       // Displacement
       displ[0]=x[0]-xi[0];
       displ[1]=x[1]-xi[1];

       // Velocity (1st deriv)
       el_pt->interpolated_dxdt(s,1,dxdt);

       some_file << x[0] << " " << x[1] << " " 
                 << displ[0] << " " << displ[1] << " "  
                 << dxdt[0] << " " << dxdt[1] << " "  
                 << std::endl;
      }
    }
  }
 some_file.close(); 

}



//========================================================================
/// Dump the solution
//========================================================================
template<class ELEMENT, class TIMESTEPPER>
void DiskShockWaveProblem<ELEMENT,TIMESTEPPER>::dump_it(ofstream& dump_file)
{
 // Call generic dump()
 Problem::dump(dump_file);
}



//========================================================================
/// Read solution from disk
//========================================================================
template<class ELEMENT, class TIMESTEPPER>
void DiskShockWaveProblem<ELEMENT,TIMESTEPPER>::restart(ifstream& restart_file)
{
 // Read generic problem data
 Problem::read(restart_file);
}



//==================================================================
/// Run the problem; specify case_number to label output directory
//==================================================================
template<class ELEMENT, class TIMESTEPPER>
void DiskShockWaveProblem<ELEMENT,TIMESTEPPER>::run(
 const unsigned& case_number)
{

 // If there's a command line argument, run the validation (i.e. do only 
 // 3 timesteps; otherwise do a few cycles
 unsigned nstep=400;
 if (CommandLineArgs::Argc!=1)
  {
   nstep=3;
  }

 // Define output directory
 char dirname[100];
 sprintf(dirname,"RESLT%i",case_number);
 Doc_info.set_directory(dirname);

 // Step number
 Doc_info.number()=0;

 // Open trace file
 char filename[100];   
 sprintf(filename,"%s/trace.dat",Doc_info.directory().c_str());
 Trace_file.open(filename);

 // Set up trace nodes as the nodes on boundary 1 (=curved boundary) in 
 // the original mesh (they exist under any refinement!) 
 unsigned nnod0=solid_mesh_pt()->nboundary_node(0);
 unsigned nnod1=solid_mesh_pt()->nboundary_node(1);
 Trace_file << "VARIABLES=\"time\"";
 for (unsigned j=0;j<nnod0;j++)
  {
   Trace_file << ", \"radial node " << j << "\" ";
  }
 for (unsigned j=0;j<nnod1;j++)
  {
   Trace_file << ", \"azimuthal node " << j << "\" ";
  }
 Trace_file << std::endl;



//  // Restart?
//  //---------

//  // Pointer to restart file
//  ifstream* restart_file_pt=0;

//  // No restart
//  //-----------
//  if (CommandLineArgs::Argc==1)
//   {
//    cout << "No restart" << std::endl;
//   }
//  // Restart
//  //--------
//  else if (CommandLineArgs::Argc==2)
//   {
//    // Open restart file
//    restart_file_pt=new ifstream(CommandLineArgs::Argv[1],ios_base::in);
//    if (restart_file_pt!=0)
//     {
//      cout << "Have opened " << CommandLineArgs::Argv[1] << 
//       " for restart. " << std::endl;
//     }
//    else
//     {
//      cout << "ERROR while trying to open " << CommandLineArgs::Argv[1] << 
//       " for restart." << std::endl;
//     }
//    // Do the actual restart
//    pause("need to do the actual restart");
//    //problem.restart(*restart_file_pt);
//   }
//  // More than one restart file specified?
//  else 
//   {
//    cout << "Can only specify one input file " << std::endl;
//    cout << "You specified the following command line arguments: " << std::endl;
//    CommandLineArgs::output();
//    //assert(false);
//   }


 // Initial parameter values
 Global_Physical_Variables::P = 0.1; 

 // Initialise time
 double time0=0.0;
 time_pt()->time()=time0;

 // Set initial timestep
 double dt=0.01; 

 // Impulsive start
 assign_initial_values_impulsive(dt); 
 
 // Doc initial state
 doc_solution();
 Doc_info.number()++;

 // First step without adaptivity
 unsteady_newton_solve(dt); 
 doc_solution();
 Doc_info.number()++;

 //Timestepping loop for subsequent steps with adaptivity
 unsigned max_adapt=1;
 for(unsigned i=1;i<nstep;i++)
  {
   unsteady_newton_solve(dt,max_adapt,false);
   doc_solution();
   Doc_info.number()++;
  }

}





//======================================================================
/// Driver for simple elastic problem
//======================================================================
int main(int argc, char* argv[])
{

 // Store command line arguments
 CommandLineArgs::setup(argc,argv);

 //Initialise physical parameters
 Global_Physical_Variables::E = 1.0; // ADJUST 
 Global_Physical_Variables::Nu = 0.3; // ADJUST

  // "Big G" Linear constitutive equations:
  Global_Physical_Variables::Constitutive_law_pt = 
   new GeneralisedHookean(&Global_Physical_Variables::Nu,
                          &Global_Physical_Variables::E);
 
  //Set up the problem:
  unsigned case_number=0;


 // Pure displacement formulation
  {
   cout << "Running case " << case_number 
        << ": Pure displacement formulation" << std::endl;
   DiskShockWaveProblem<RefineableQPVDElement<2,3>, Newmark<1> > problem;
   problem.run(case_number);
   case_number++;
  }
   
 // Pressure-displacement with Crouzeix Raviart-type pressure
  {
   cout << "Running case " << case_number 
        << ": Pressure/displacement with Crouzeix-Raviart pressure" << std::endl;
   DiskShockWaveProblem<RefineableQPVDElementWithPressure<2>, Newmark<1> >
    problem;
   problem.run(case_number);
   case_number++;
  }


  // Pressure-displacement with Taylor-Hood-type pressure
  {
   cout << "Running case " << case_number 
        << ": Pressure/displacement with Taylor-Hood pressure" << std::endl;
   DiskShockWaveProblem<RefineableQPVDElementWithContinuousPressure<2>, 
    Newmark<1> > problem;
   problem.run(case_number);
   case_number++;
  }

 
 // Clean up 
 delete Global_Physical_Variables::Constitutive_law_pt;
 Global_Physical_Variables::Constitutive_law_pt=0;

}

---