Qelements.h

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 1.0; svn revision $LastChangedRevision$
//LIC//
//LIC// $LastChangedDate$
//LIC// 
//LIC// Copyright (C) 2006-2016 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//====================================================================
//Header functions for classes that define Qelements
//Include guards to prevent multiple inclusion of the header
#ifndef OOMPH_QELEMENT_HEADER
#define OOMPH_QELEMENT_HEADER

// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
  #include <oomph-lib-config.h>
#endif


#ifdef OOMPH_HAS_MPI
#include "mpi.h"
#endif

//oomph-lib headers
#include "Vector.h"
#include "shape.h"
#include "integral.h"
#include "timesteppers.h"
#include "elements.h"
#include "macro_element.h"

#include "Qelement_face_coordinate_translation_schemes.h"



namespace oomph
{


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

//========================================================================
/// Empty base class for Qelements (created so that 
/// we can use dynamic_cast<>() to figure out if a an element
/// is a Qelement (from a purely geometric point of view).
//========================================================================
 class QElementGeometricBase : public virtual FiniteElement
{

  public:

 /// Empty default constructor
 QElementGeometricBase(){} 

 /// Broken copy constructor
 QElementGeometricBase(const QElementGeometricBase&) 
  { 
   BrokenCopy::broken_copy("QElementGeometricBase");
  } 
 
 /// Broken assignment operator
//Commented out broken assignment operator because this can lead to a conflict warning
//when used in the virtual inheritence hierarchy. Essentially the compiler doesn't
//realise that two separate implementations of the broken function are the same and so,
//quite rightly, it shouts.
 /*void operator=(const QElementGeometricBase&) 
  {
   BrokenCopy::broken_assign("QElementGeometricBase");
   }*/


};


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


//========================================================================
/// Base class for Qelements 
//========================================================================
class QElementBase : public virtual QElementGeometricBase
{
  public:

 /// Constructor: Initialise pointers to macro element reference coords
 QElementBase() :  S_macro_ll_pt(0), S_macro_ur_pt(0) 
  {}

 /// Broken copy constructor
 QElementBase(const QElementBase&) 
  { 
   BrokenCopy::broken_copy("QElementBase");
  } 
 
 /// Broken assignment operator
 /*void operator=(const QElementBase&) 
  {
   BrokenCopy::broken_assign("QElementBase");
   }*/

 /// Destructor: Kill storage for macro element reference coords
 virtual ~QElementBase()
  {
   // Can be deleted blindly as they were nulled initially
   delete S_macro_ll_pt;
   S_macro_ll_pt=0;
   delete S_macro_ur_pt;
   S_macro_ur_pt=0;
  }

 ///Check whether the local coordinate are valid or not 
 bool local_coord_is_valid(const Vector<double> &s) 
  {
   unsigned ncoord = dim();
   for(unsigned i=0;i<ncoord;i++)
    {
     // We're outside 
     if((s[i] - s_max() >  0.0) ||
        (s_min() - s[i] >  0.0))
      {
       return false;
      }
    }
   return true;
  }
 
 /// \short Adjust local coordinates so that they're located inside
 /// the element
 void move_local_coord_back_into_element(Vector<double> &s) const
 {
  unsigned ncoord = dim();
  for(unsigned i=0;i<ncoord;i++)
   {
    // Adjust to move it onto the boundary
    if (s[i] > s_max() ) s[i] = s_max();
    if (s[i] < s_min() ) s[i] = s_min();
   }
 }
 

 /// \short Set pointer to macro element also sets up storage for the
 /// reference coordinates and initialises them
 virtual void set_macro_elem_pt(MacroElement* macro_elem_pt)
  {
   // Get the spatial dimension (= number of local and macro element coords)
   unsigned n_dim=dim();

   // Create storage if none has been allocated
   if(S_macro_ll_pt==0)
    {
     S_macro_ll_pt=new Vector<double>(n_dim);
    }
   //Otherwise resize the allocated storage
   else
    {
     S_macro_ll_pt->resize(n_dim);
    }
   
   //Create storage if none has been allocated
   if(S_macro_ur_pt==0)
    {
     S_macro_ur_pt=new Vector<double>(n_dim);
    }
   //Otherwise resize the allocated storage
   else
    {
     S_macro_ur_pt->resize(n_dim);
    }

   // Initialise the vertex coordinates in the macro element 
   // Default: The element is unrefined and hence its vertices are those
   // of the macro element itself
   for (unsigned i=0;i<n_dim;i++)
    {
     s_macro_ll(i)=-1.0; 
     s_macro_ur(i)= 1.0;
    }

   /// Call the corresponding function in the FiniteElement base class
   FiniteElement::set_macro_elem_pt(macro_elem_pt);
  }
   

 /// \short Access fct to the i-th coordinate of the element's
 /// "lower left" vertex in the associated MacroElement
 double& s_macro_ll(const unsigned& i)
  {
#ifdef PARANOID
   if (S_macro_ll_pt==0)
    {
     throw OomphLibError("S_macro_ll_pt has not been set\n",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif   
   return (*S_macro_ll_pt)[i];
  }



 /// \short Access fct to the i-th coordinate of the element's
 /// "upper right" vertex in the associated MacroElement
 double& s_macro_ur(const unsigned& i)
  {
#ifdef PARANOID
   if (S_macro_ur_pt==0)
    {
     throw OomphLibError("S_macro_ur_pt has not been set\n",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif   
   return (*S_macro_ur_pt)[i];
  }



 /// \short Access fct to the i-th coordinate of the element's
 /// "lower left" vertex in the associated MacroElement. (const version)
 double s_macro_ll(const unsigned& i) const
  {
#ifdef PARANOID
   if (S_macro_ll_pt==0)
    {
     throw OomphLibError("S_macro_ll_pt has not been set\n",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif   
   return (*S_macro_ll_pt)[i];
  }



 /// \short Access fct to the i-th coordinate of the element's
 /// "upper right" vertex in the associated MacroElement. (const version)
 double s_macro_ur(const unsigned& i) const
  {
#ifdef PARANOID
   if (S_macro_ur_pt==0)
    {
     throw OomphLibError("S_macro_ur_pt has not been set\n",
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
#endif   
   return (*S_macro_ur_pt)[i];
  }

 /// \short Global coordinates as function of local coordinates.
 /// using the macro element representation
 void get_x_from_macro_element(const Vector<double> &s, 
                               Vector<double> &x) const
  {
   //Check that there is a macro element
   if(Macro_elem_pt==0) 
    {
     throw OomphLibError(
      "Macro Element pointer not set in this element\n",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }
   
   //Use macro element representation
   unsigned el_dim = dim();
   Vector<double> s_macro(el_dim);
   for(unsigned i=0;i<el_dim;i++)
    {
     s_macro[i]=s_macro_ll(i)+0.5*(s[i]+1.0)*(s_macro_ur(i)-s_macro_ll(i));
    }
   Macro_elem_pt->macro_map(s_macro,x);
  }
 
 /// \short Global coordinates as function of local coordinates
 /// at previous time "level" t (t=0: present; t>0: previous)
 /// using the macro element representation
 void get_x_from_macro_element(const unsigned& t, const Vector<double>& s, 
                               Vector<double>& x) 
  {
   //Check that there is a macro element
   if(Macro_elem_pt==0) 
    {
     throw OomphLibError(
      "Macro Element pointer not set in this element\n",
      OOMPH_CURRENT_FUNCTION,
      OOMPH_EXCEPTION_LOCATION);
    }

   //Use the macro element representation
   unsigned el_dim = dim();
   Vector<double> s_macro(el_dim);
   for(unsigned i=0;i<el_dim;i++)
    {
     s_macro[i]=s_macro_ll(i)+0.5*(s[i]+1.0)*(s_macro_ur(i)-s_macro_ll(i));
    }
   Macro_elem_pt->macro_map(t,s_macro,x);
  }

 /// Return number of nodes on one face of the element. Always
 /// nnode_1d^(el_dim - 1).
 unsigned nnode_on_face() const
  {
   // c++ doesn't have pow(int, int) so we have to use all these casts...
   return static_cast<unsigned>(std::pow(static_cast<double>(nnode_1d()), 
                                          static_cast<double>(dim()-1)));
  }

 /// It's a Q element!
 ElementGeometry::ElementGeometry element_geometry() const
 {
  return ElementGeometry::Q;
 }

  private:

 /// Pointer to vector of lower left vertex coords. in macro element
 Vector<double>* S_macro_ll_pt;
 
 /// Pointer to vector of upper right vertex coords. in macro element
 Vector<double>* S_macro_ur_pt;

};


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



//========================================================================
/// Base class for Solid Qelements 
//========================================================================
class QSolidElementBase : public virtual QElementBase,
                          public virtual SolidFiniteElement
{
  public:

 /// Constructor: Empty
 QSolidElementBase(){};

 /// Broken copy constructor
 QSolidElementBase(const QSolidElementBase&) 
  { 
   BrokenCopy::broken_copy("QSolidElementBase");
  } 
 
 /// Broken assignment operator
 /*void operator=(const QSolidElementBase&) 
  {
   BrokenCopy::broken_assign("QSolidElementBase");
   }*/

 /// \short Set pointer to MacroElement -- overloads generic version
 /// in RefineableQElement<2> and uses the MacroElement 
 /// also as the default for the "undeformed" configuration.
 /// This assignment can/must be overwritten with 
 /// set_undeformed_macro_elem_pt(...) if the deformation of 
 /// the solid body is driven by a deformation of the
 /// "current" Domain/MacroElement representation of it's boundary.
 virtual void set_macro_elem_pt(MacroElement* macro_elem_pt)
  {
   // Call the general Q version which sets up the storage
   // for the reference coordinates
   QElementBase::set_macro_elem_pt(macro_elem_pt);
   // Store pointer to macro element that represents the exact
   // undeformed geomtry
   set_undeformed_macro_elem_pt(macro_elem_pt);
  }
 
 /// \short Set pointers to "current" and "undeformed" MacroElements.
 virtual void set_macro_elem_pt(MacroElement* macro_elem_pt,
                                MacroElement* undeformed_macro_elem_pt)
  {
   // Call the general Q version which sets up the storage
   // for the reference coordinates
   QElementBase::set_macro_elem_pt(macro_elem_pt);
   // Store pointer to macro element that represents the exact
   // undeformed geomtry
   set_undeformed_macro_elem_pt(undeformed_macro_elem_pt);
  }

 /// \short Eulerian and Lagrangian coordinates as function of the
 /// local coordinates: The Eulerian position is returned in 
 /// FE-interpolated form (\c x_fe) and then in the form obtained
 /// from the "current" MacroElement representation (if it exists -- if not,
 /// \c x is the same as \c x_fe). This allows the Domain/MacroElement-
 /// based representation to be used to apply displacement boundary
 /// conditions exactly. Ditto for the Lagrangian coordinates returned
 /// in xi_fe and xi.
void get_x_and_xi(const Vector<double>& s, 
                  Vector<double>& x_fe,
                  Vector<double>& x,
                  Vector<double>& xi_fe,
                  Vector<double>& xi) const
 {
  // Lagrangian coordinate: Directly from 
  // underlying FE representation
  unsigned n_xi = xi_fe.size();
  for(unsigned i=0;i<n_xi;i++) {xi_fe[i] = interpolated_xi(s,i);}

  // Lagrangian coordinate from FE representation again
  if(Undeformed_macro_elem_pt==0)
   {
    unsigned n_xi = xi.size();
    for(unsigned i=0;i<n_xi;i++) {xi[i] = xi_fe[i];}
   }
  // ...or refer to the "undeformed" MacroElement if it exists.
  else 
   {
    unsigned el_dim = dim();
    Vector<double> s_macro(el_dim);
    for(unsigned i=0;i<el_dim;i++)
     {
      s_macro[i]=s_macro_ll(i)+0.5*(s[i]+1.0)*(s_macro_ur(i)-s_macro_ll(i));
     }
    Undeformed_macro_elem_pt->macro_map(s_macro,xi);
   }
  
  
  // Eulerian coordinate directly from  underlying FE representation
  unsigned n_x=x_fe.size();
  for(unsigned i=0;i<n_x;i++) {x_fe[i] = interpolated_x(s,i);}
  
  // Eulerian coordinate from FE representation again:
  if(Macro_elem_pt==0)
   {
    for(unsigned i=0;i<n_x;i++) {x[i] = x_fe[i];}
   }
  // or refer to the "current" MacroElement if it exists.
  else 
   {
    unsigned el_dim = dim();
    Vector<double> s_macro(el_dim);
    for(unsigned i=0;i<el_dim;i++)
     {
      s_macro[i]=s_macro_ll(i)+0.5*(s[i]+1.0)*(s_macro_ur(i)-s_macro_ll(i));
     }
    Macro_elem_pt->macro_map(s_macro,x);
   }
 }

};




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



//=======================================================================
/// General QElement class
///
/// Empty, just establishes the template parameters
//=======================================================================
template<unsigned DIM, unsigned NNODE_1D> 
 class QElement 
{
};

//=======================================================================
/// Base class for all line elements 
//=======================================================================
class LineElementBase : public virtual QElementBase
{
  public:

 /// Constructor. Empty
 LineElementBase() {}

 /// \short Number of vertex nodes in the element
 virtual unsigned nvertex_node() const=0;

 /// \short Pointer to the j-th vertex node in the element
 virtual Node* vertex_node_pt(const unsigned& j) const=0;

};

//=======================================================================
/// General QElement class specialised to one spatial dimension
//=======================================================================
template<unsigned NNODE_1D>
class QElement<1,NNODE_1D> : public virtual LineElementBase
{
  private:

 /// \short Default integration rule: Gaussian integration of same 'order' 
 /// as the element
 //This is sort of optimal, because it means that the integration is exact
 //for the shape functions. Can overwrite this in specific element defintion.
 static Gauss<1,NNODE_1D> Default_integration_scheme;

public:

 /// Constructor
 QElement() 
  {
   //There are NNODE_1D nodes in this element
   this->set_n_node(NNODE_1D);
   //Set the dimensions of the element and the nodes, by default, both 1D
   this->set_dimension(1);
   //Assign pointer to default (full) integration_scheme
   this->set_integration_scheme(&Default_integration_scheme);
  }

 /// Broken copy constructor
 QElement(const QElement&) 
  { 
   BrokenCopy::broken_copy("QElement");
  } 
 
 /// Broken assignment operator
 /*void operator=(const QElement&) 
  {
   BrokenCopy::broken_assign("QElement");
   }*/

 /// Calculate the geometric shape functions at local coordinate s
 void shape(const Vector<double> &s, Shape &psi) const;

 /// \short Compute the  geometric shape functions and 
 /// derivatives w.r.t. local coordinates at local coordinate s
 void dshape_local(const Vector<double> &s, Shape &psi, DShape &dpsids) 
  const;

 /// \short Compute the geometric shape functions, derivatives and
 /// second derivatives w.r.t. local coordinates at local coordinate s 
 /// d2psids(i,0) = \f$ d^2 \psi_j / d s^2 \f$
 void d2shape_local(const Vector<double> &s, Shape &psi, DShape &dpsids,
                    DShape &d2psids) const;

 /// \short Overload the template-free interface for the calculation of
 /// inverse jacobian matrix. This is a one-dimensional element, so
 /// use the 1D version.
 double invert_jacobian_mapping(const DenseMatrix<double> &jacobian,
                                DenseMatrix<double> &inverse_jacobian) const
  {return FiniteElement::invert_jacobian<1>(jacobian,inverse_jacobian);}

 /// Min. value of local coordinate
 double s_min() const {return -1.0;}

 /// Max. value of local coordinate
 double s_max() const {return  1.0;}

 /// \short Number of vertex nodes in the element
 unsigned nvertex_node() const
  {return 2;}
 
 /// \short Pointer to the j-th vertex node in the element
 Node* vertex_node_pt(const unsigned& j) const
  {
   unsigned n_node_1d=nnode_1d();
   Node* nod_pt;
   switch(j)
    {
    case 0:
     nod_pt=node_pt(0);
     break;
    case 1:
     nod_pt=node_pt(n_node_1d-1);
     break;
    default:
     std::ostringstream error_message;
     error_message << "Vertex node number is " << j << 
      " but must be from 0 to 1\n";
     
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   return nod_pt;
  }
 
 /// Get local coordinates of node j in the element; vector sets its own size
 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
  {
   s.resize(1);
   s[0]=this->s_min()+double(j)/double(NNODE_1D-1)*
    (this->s_max()-this->s_min());
  }

 
 /// Get the local fraction of node j in the element
 void local_fraction_of_node(const unsigned &j,
                             Vector<double> &s_fraction)
  {
   s_fraction.resize(1);
   s_fraction[0] = double(j)/double(NNODE_1D-1);
  }


 /// This function returns the local fraction of all nodes at the n-th 
 /// position in a one dimensional expansion along the i-th local coordinate
 inline double local_one_d_fraction_of_node(
  const unsigned &n1d,const unsigned &i)
  {
   //It's just the value of the node divided by the number of 1-D nodes
   return double(n1d)/double(NNODE_1D-1);
  }

 /// Get the node at the specified local coordinate
 Node* get_node_at_local_coordinate(const Vector<double> &s) const;

 /// Number of nodes along each element edge
 unsigned nnode_1d() const {return NNODE_1D;}

 /// \short Return the number of actual plot points for paraview
 /// plot with parameter nplot.
 unsigned nplot_points_paraview(const unsigned& nplot) const
  {
   return nplot;
  }
 
 /// \short Return the number of local sub-elements for paraview plot with 
 /// parameter nplot.
 unsigned nsub_elements_paraview(const unsigned& nplot) const 
  {
   return (nplot-1);   
  }

 /// \short Fill in the offset information for paraview plot.
 /// Needs to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 void write_paraview_output_offset_information(std::ofstream& file_out,
                                               const unsigned& nplot,
                                               unsigned& counter) const
  {
   // Number of local elements we want to plot over
   unsigned plot=nsub_elements_paraview(nplot);
   
   // loops over the i-th local element in parent element
   for(unsigned i=0;i<plot;i++)
    {
     file_out << i+counter << " "
              << i+1+counter 
              << std::endl;
    }
   counter+=nplot_points_paraview(nplot);
  }

 /// \short Return the paraview element type.
 /// Needs to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 /// Use type "VTK_LINE" (== 3) for 2D quad elements
 void write_paraview_type(std::ofstream& file_out,
                          const unsigned& nplot) const
  {
   unsigned local_loop=nsub_elements_paraview(nplot);
   for(unsigned i=0;i<local_loop;i++)
    {
     file_out << "3" << std::endl;
    }
  }
 
 /// \short Return the offsets for the paraview sub-elements. Needs 
 /// to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 void write_paraview_offsets(std::ofstream& file_out,
                             const unsigned& nplot, 
                             unsigned& offset_sum) const
  {
   // Loop over all local elements and add its offset to the overall offset_sum
   unsigned local_loop=nsub_elements_paraview(nplot);
   for(unsigned i=0;i<local_loop;i++)
    {
     offset_sum+=2;
     file_out << offset_sum << std::endl;
    }
  }

 /// Output 
 void output(std::ostream &outfile);

 /// Output at n_plot points
 void output(std::ostream &outfile, const unsigned &n_plot);

 /// C-style output 
 void output(FILE* file_pt);

 /// C_style output at n_plot points
 void output(FILE* file_pt, const unsigned &n_plot);


 /// \short  Get cector of local coordinates of plot point i (when plotting 
 /// nplot points in each "coordinate direction).
 void get_s_plot(
  const unsigned& i, 
  const unsigned& nplot,
  Vector<double>& s,
  const bool& use_equally_spaced_interior_sample_points=false) const
  {
   if (nplot>1)
    {
     s[0]=-1.0+2.0*double(i)/double(nplot-1);     
    if (use_equally_spaced_interior_sample_points)
     {
      double range=2.0;
      double dx_new=range/double(nplot);
      double range_new=double(nplot-1)*dx_new;
      s[0]=-1.0+0.5*dx_new+range_new*(1.0+s[0])/range;
     }
    }
   else
    {
     s[0]=0.0;
    }
  }

 /// \short Return string for tecplot zone header (when plotting 
 /// nplot points in each "coordinate direction)
 std::string tecplot_zone_string(const unsigned& nplot) const
 {
  std::ostringstream header;
  header << "ZONE I=" << nplot << "\n";
  return header.str();
 }

 /// \short Return total number of plot points (when plotting 
 /// nplot points in each "coordinate direction)
 unsigned nplot_points(const unsigned& nplot) const
  {
   unsigned DIM=1;
   unsigned np=1;
   for (unsigned i=0;i<DIM;i++) {np*=nplot;}
   return np;
  }

 /// Get the number of the ith node on face face_index in the bulk node
 /// vector. 
 unsigned get_bulk_node_number(const int& face_index,
                               const unsigned& i) const
  {
   face_node_number_error_check(i);

   if(face_index == -1) {return 0;}
   else if(face_index == +1) {return nnode_1d() - 1;}
   else
    {
     std::string err = "Face index should be in {-1, +1}.";
     throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                         OOMPH_CURRENT_FUNCTION);
    }
  }

 /// Get the sign of the outer unit normal on the face given by face_index.
 int face_outer_unit_normal_sign(const int& face_index) const
  {
#ifdef PARANOID
   if(std::abs(face_index) != 1)
    {
     std::string err = "Face index should be in {-1, +1}.";
     throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                         OOMPH_CURRENT_FUNCTION);
    }
#endif
   return face_index;
  }

 /// Get a pointer to the function mapping face coordinates to bulk coordinates
 CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt
 (const int& face_index) const
 {
  if(face_index == 1) {return &QElement1FaceToBulkCoordinates::face1;}
  else if(face_index == -1) {return &QElement1FaceToBulkCoordinates::face0;}
  else
   {
    std::string err = "Face index should be in {-1, +1}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   } 
 }

 /// Get a pointer to the derivative of the mapping from face to bulk
 /// coordinates.
 BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt
 (const int& face_index) const
 {
  if(face_index == 1) {return &QElement1BulkCoordinateDerivatives::faces0;}
  else if (face_index == -1) {return &QElement1BulkCoordinateDerivatives::faces0;}
  else
   {
    std::string err = "Face index should be in {-1, +1}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   } 
 }

};

//=======================================================================
/// Base class for all quad elements 
//=======================================================================
class QuadElementBase : public virtual QElementBase
{
  public:

 /// Constructor. Empty
 QuadElementBase() {}

 /// \short Number of vertex nodes in the element
 virtual unsigned nvertex_node() const=0;

 /// \short Pointer to the j-th vertex node in the element
 virtual Node* vertex_node_pt(const unsigned& j) const=0;

};

//=======================================================================
/// General QElement class specialised to two spatial dimensions
//=======================================================================
template<unsigned NNODE_1D>
class QElement<2,NNODE_1D> : public virtual QuadElementBase
{
  private:
  
 /// \short Default integration rule: Gaussian integration of same 'order' 
 /// as the element
 //N.B. This is sort of optimal, because it means that the integration is exact
 //for the shape functions. Can overwrite this in specific element defintion 
 static Gauss<2,NNODE_1D> Default_integration_scheme; 
 
public: 

 ///Constructor
 QElement() 
  {
   //There are NNODE_1D*NNODE_1D nodes in this element
   this->set_n_node(NNODE_1D*NNODE_1D);
   //Set the dimensions of the element and the nodes, by default, both 2D
   set_dimension(2);
   //Assign default (full) spatial integration scheme
   set_integration_scheme(&Default_integration_scheme);   
  }
  
 /// Broken copy constructor
 QElement(const QElement&) 
  { 
   BrokenCopy::broken_copy("QElement");
  } 
 
 /// Broken assignment operator
 /*void operator=(const QElement&) 
  {
   BrokenCopy::broken_assign("QElement");
   }*/

 /// Calculate the geometric shape functions at local coordinate s
 void shape(const Vector<double> &s, Shape &psi) const;

 /// \short Compute the  geometric shape functions and 
 /// derivatives w.r.t. local coordinates at local coordinate s
 void dshape_local(const Vector<double> &s, Shape &psi, DShape &dpsids) 
  const;

 /// \short Compute the geometric shape functions, derivatives and
 /// second derivatives w.r.t. local coordinates at local coordinate s 
 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$ 
 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$ 
 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$ 
 void d2shape_local(const Vector<double> &s, Shape &psi, DShape &dpsids,
                    DShape &d2psids) const;

 /// \short Overload the template-free interface for the calculation of
 /// inverse jacobian matrix. This is a two-dimensional element, so use
 /// the two-d version.
 double invert_jacobian_mapping(const DenseMatrix<double> &jacobian,
                                DenseMatrix<double> &inverse_jacobian) const
  {return FiniteElement::invert_jacobian<2>(jacobian,inverse_jacobian);}
 
 /// Min. value of local coordinate
 double s_min() const {return -1.0;}

 /// Max. value of local coordinate
 double s_max() const {return  1.0;}


 /// \short Number of vertex nodes in the element
 unsigned nvertex_node() const
  {return 4;}

 /// \short Pointer to the j-th vertex node in the element
 Node* vertex_node_pt(const unsigned& j) const
  {
   unsigned n_node_1d=nnode_1d();
   Node* nod_pt;
   switch(j)
    {
    case 0:
     nod_pt=node_pt(0);
     break;
    case 1:
     nod_pt=node_pt(n_node_1d-1);
     break;
    case 2:
     nod_pt=node_pt(n_node_1d*(n_node_1d-1));
     break;
    case 3:
     nod_pt=node_pt(n_node_1d*n_node_1d-1);
     break;
    default:
     std::ostringstream error_message;
     error_message << "Vertex node number is " << j << 
      " but must be from 0 to 3\n";
     
     throw OomphLibError(error_message.str(),
                         OOMPH_CURRENT_FUNCTION,
                         OOMPH_EXCEPTION_LOCATION);
    }
   return nod_pt;
  }
 

 /// Get local coordinates of node j in the element; vector sets its own size
 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
  {
   s.resize(2);
   unsigned j0=j%NNODE_1D;
   unsigned j1=j/NNODE_1D;
   const double S_min = this->s_min();
   const double S_range = this->s_max() - S_min;
   s[0]=S_min+double(j0)/double(NNODE_1D-1)*S_range;
   s[1]=S_min+double(j1)/double(NNODE_1D-1)*S_range;
  }

 /// Get the local fraction of node j in the element
 void local_fraction_of_node(const unsigned &j,
                             Vector<double> &s_fraction)
  {
   s_fraction.resize(2);
   unsigned j0 = j%NNODE_1D;
   unsigned j1 = j/NNODE_1D;
   s_fraction[0] = double(j0)/double(NNODE_1D-1);
   s_fraction[1] = double(j1)/double(NNODE_1D-1);
  }

 /// This function returns the local fraction of ant nodes in the n-th positoin
 /// in a one dimensional expansion along the i-th local coordinate
 inline double local_one_d_fraction_of_node(const unsigned &n1d,
                                            const unsigned &i)
  {
   //It's just the value of the node divided by the number of 1-D nodes
   return double(n1d)/double(NNODE_1D-1);
  }

 /// Get the node at the specified local coordinate
 Node* get_node_at_local_coordinate(const Vector<double> &s) const;
 
 /// Number of nodes along each element edge
 unsigned nnode_1d() const {return NNODE_1D;}

 /// \short Return the number of actual plot points for paraview
 /// plot with parameter nplot.
 unsigned nplot_points_paraview(const unsigned& nplot) const
  {
   return nplot*nplot;
  }
 
 /// \short Return the number of local sub-elements for paraview plot with 
 /// parameter nplot.
 unsigned nsub_elements_paraview(const unsigned& nplot) const 
  {
   return (nplot-1)*(nplot-1);   
  }

 /// \short Fill in the offset information for paraview plot.
 /// Needs to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 void write_paraview_output_offset_information(std::ofstream& file_out,
                                               const unsigned& nplot,
                                               unsigned& counter) const
  {
   // Number of local elements we want to plot over
   unsigned plot=nsub_elements_paraview(nplot);
   
   // loops over the i-th local element in parent element
   for(unsigned i=0;i<plot;i++)
    {
     unsigned d=(i-(i%(nplot-1)))/(nplot-1);
     
     file_out << i%(nplot-1)+d*nplot+counter << " "
              << i%(nplot-1)+1+d*nplot+counter << " "
              << i%(nplot-1)+1+(d+1)*nplot+counter << " "
              << i%(nplot-1)+(d+1)*nplot+counter
              << std::endl;
    }
   counter+=nplot_points_paraview(nplot);
  }

 /// \short Return the paraview element type.
 /// Needs to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 /// Use type "VTK_QUAD" (== 9) for 2D quad elements
 void write_paraview_type(std::ofstream& file_out,
                          const unsigned& nplot) const
  {
   unsigned local_loop=nsub_elements_paraview(nplot);
   for(unsigned i=0;i<local_loop;i++)
    {
     file_out << "9" << std::endl;
    }
  }
 
 /// \short Return the offsets for the paraview sub-elements. Needs 
 /// to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 void write_paraview_offsets(std::ofstream& file_out,
                             const unsigned& nplot, 
                             unsigned& offset_sum) const
  {
   // Loop over all local elements and add its offset to the overall offset_sum
   unsigned local_loop=nsub_elements_paraview(nplot);
   for(unsigned i=0;i<local_loop;i++)
    {
     offset_sum+=4;
     file_out << offset_sum << std::endl;
    }
  }

 /// Output
 void output(std::ostream &outfile);

 /// Output at n_plot points
 void output(std::ostream &outfile, const unsigned &n_plot);

 /// C-style output 
 void output(FILE* file_pt);

 /// C_style output at n_plot points
 void output(FILE* file_pt, const unsigned &n_plot);


 /// \short  Get cector of local coordinates of plot point i (when plotting 
 /// nplot points in each "coordinate direction).
 void get_s_plot(
  const unsigned& i, 
  const unsigned& nplot,
  Vector<double>& s,
  const bool& use_equally_spaced_interior_sample_points=false) const
  {
   if (nplot>1)
    {
     unsigned i0=i%nplot;
     unsigned i1=(i-i0)/nplot;
     
     s[0]=-1.0+2.0*double(i0)/double(nplot-1);
     s[1]=-1.0+2.0*double(i1)/double(nplot-1); 
    if (use_equally_spaced_interior_sample_points)
     {
      double range=2.0;
      double dx_new=range/double(nplot);
      double range_new=double(nplot-1)*dx_new;
      s[0]=-1.0+0.5*dx_new+range_new*(1.0+s[0])/range;
      s[1]=-1.0+0.5*dx_new+range_new*(1.0+s[1])/range;
     }
    }
   else
    {
     s[0]=0.0;
     s[1]=0.0;
    }
  }

 /// \short Return string for tecplot zone header (when plotting 
 /// nplot points in each "coordinate direction)
 std::string tecplot_zone_string(const unsigned& nplot) const
  {
   std::ostringstream header;
   header << "ZONE I=" << nplot << ", J=" << nplot << "\n";
   return header.str();
  }

 /// Return total number of plot points (when plotting 
 /// nplot points in each "coordinate direction)
 unsigned nplot_points(const unsigned& nplot) const
  {
   unsigned DIM=2;
   unsigned np=1;
   for (unsigned i=0;i<DIM;i++) {np*=nplot;}
   return np;
  };

 /// Get the number of the ith node on face face_index (in the bulk node
 /// vector). 
 unsigned get_bulk_node_number(const int& face_index,
                               const unsigned& i) const
 {
  face_node_number_error_check(i);

  const unsigned nn1d = nnode_1d();

  if(face_index == -1) {return i*nn1d;}
  else if(face_index == +1) {return nn1d*i + nn1d-1;}
  else if(face_index == -2) {return i;}
  else if(face_index == +2) {return nn1d*(nn1d-1) + i;}
  else
   {
    std::string err = "Face index should be in {-1, +1, -2, +2}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   }
 }

 /// Get the sign of the outer unit normal on the face given by face_index.
 int face_outer_unit_normal_sign(const int& face_index) const
 {
  if(face_index < 0) {return 1;}
  else if(face_index > 0) {return -1;}
  else
   {
    std::string err = "Face index should be one of {-1, +1, -2, +2}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   }
 }

 /// Get a pointer to the function mapping face coordinates to bulk coordinates
 CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt
 (const int& face_index) const
 {
  if(face_index == 1) {return &QElement2FaceToBulkCoordinates::face2;}
  else if (face_index == -1) {return &QElement2FaceToBulkCoordinates::face0;}
  else if (face_index == -2) {return &QElement2FaceToBulkCoordinates::face1;}
  else if (face_index == 2) {return &QElement2FaceToBulkCoordinates::face3;}
  else
   {
    std::string err = "Face index should be in {-1, +1}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   } 
 }

 /// Get a pointer to the derivative of the mapping from face to bulk
 /// coordinates.
 BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt
 (const int& face_index) const
 {
  if(face_index == 1) {return &QElement2BulkCoordinateDerivatives::faces0;}
  else if (face_index == -1) {return &QElement2BulkCoordinateDerivatives::faces0;}
  else if (face_index == -2) {return &QElement2BulkCoordinateDerivatives::faces1;}
  else if (face_index == 2) {return &QElement2BulkCoordinateDerivatives::faces1;}
  else
   {
    std::string err = "Face index should be in {-1, +1}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   } 
 }



};

//=======================================================================
/// Base class for all brick elements
//=======================================================================
class BrickElementBase : public virtual QElementBase
{
  public:

 /// Constructor. Empty
 BrickElementBase() {}

 /// \short Number of vertex nodes in the element
 virtual unsigned nvertex_node() const=0;

 /// \short Pointer to the j-th vertex node in the element
 virtual Node* vertex_node_pt(const unsigned& j) const=0;

};

//=======================================================================
///General QElement class specialised to three spatial dimensions
//=======================================================================
template<unsigned NNODE_1D>
class QElement<3,NNODE_1D> : public virtual BrickElementBase
{
  private:
  
 /// \short Default integration rule: Gaussian integration of same 'order' 
 /// as the element
 //N.B. This is sort of optimal, because it means that the integration is exact
 //for the shape functions. Can overwrite this in specific element defintion 
 static Gauss<3,NNODE_1D> Default_integration_scheme; 
 
public: 

 /// Constructor
 QElement() 
  {
   //There are NNODE_1D^3 nodes in this element
   this->set_n_node(NNODE_1D*NNODE_1D*NNODE_1D);
   //Set the dimensions of the element and the nodes, by default, both 3D
   set_dimension(3);
   //Assign default (full_ spatial integration_scheme
   set_integration_scheme(&Default_integration_scheme);
  }


 /// Broken copy constructor
 QElement(const QElement&) 
  { 
   BrokenCopy::broken_copy("QElement");
  } 
 
 /// Broken assignment operator
 /*void operator=(const QElement&) 
  {
   BrokenCopy::broken_assign("QElement");
   }*/
 
 /// Calculate the geometric shape functions at local coordinate s
 void shape(const Vector<double> &s, Shape &psi) const;

 /// \short Compute the  geometric shape functions and 
 /// derivatives w.r.t. local coordinates at local coordinate s
 void dshape_local(const Vector<double> &s, Shape &psi, DShape &dpsids) 
  const;

 /// \short Compute the geometric shape functions, derivatives and
 /// second derivatives w.r.t. local coordinates at local coordinate s. 
 /// d2psids(i,0) = \f$ \partial^2 \psi_j / \partial s_0^2 \f$ 
 /// d2psids(i,1) = \f$ \partial^2 \psi_j / \partial s_1^2 \f$ 
 /// d2psids(i,2) = \f$ \partial^2 \psi_j / \partial s_2^2 \f$ 
 /// d2psids(i,3) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_1 \f$ 
 /// d2psids(i,4) = \f$ \partial^2 \psi_j / \partial s_0 \partial s_2 \f$ 
 /// d2psids(i,5) = \f$ \partial^2 \psi_j / \partial s_1 \partial s_2 \f$
 void d2shape_local(const Vector<double> &s, Shape &psi, DShape &dpsids,
                    DShape &d2psids) const;


 /// \short Overload the template-free interface for the calculation of
 /// the inverse jacobian mapping. This is a three-dimensional element,
 /// so use the 3d version
 double invert_jacobian_mapping(const DenseMatrix<double> &jacobian,
                                DenseMatrix<double> &inverse_jacobian) const
  {return FiniteElement::invert_jacobian<3>(jacobian,inverse_jacobian);}

 /// Min. value of local coordinate
 double s_min() const {return -1.0;}

 /// Max. value of local coordinate
 double s_max() const {return  1.0;}
 
 /// \short Number of vertex nodes in the element
 unsigned nvertex_node() const
  {return 8;}

 /// \short Pointer to the j-th vertex node in the element
 Node* vertex_node_pt(const unsigned& j) const
 {
  unsigned N=nnode_1d();
  Node* nod_pt;
  switch(j)
   {
   case 0:
    nod_pt=node_pt(0);
    break;
   case 1:
    nod_pt=node_pt(N-1);
    break;
   case 2:
    nod_pt=node_pt(N*(N-1));
    break;
   case 3:
    nod_pt=node_pt(N*N-1);
    break;
   case 4:
    nod_pt=node_pt(N*N*(N-1));
    break;
   case 5:
    nod_pt=node_pt(N*N*(N-1)+(N-1));
    break;
   case 6:
    nod_pt=node_pt(N*N*N-N);
    break;
   case 7:
    nod_pt=node_pt(N*N*N-1);
    break;
   default:
    std::ostringstream error_message;
    error_message << "Vertex node number is " << j << 
     " but must be from 0 to 7\n";
    
    throw OomphLibError(error_message.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
   }
  return nod_pt;
 }

 /// Get local coordinates of node j in the element; vector sets its own size
 void local_coordinate_of_node(const unsigned& j, Vector<double>& s) const
  {
   s.resize(3);
   unsigned j0=j%NNODE_1D;
   unsigned j1=(j/NNODE_1D)%NNODE_1D;
   unsigned j2=j/(NNODE_1D*NNODE_1D);
   const double S_min = this->s_min();
   const double S_range = this->s_max() - S_min;

   s[0]=S_min+double(j0)/double(NNODE_1D-1)*S_range;
   s[1]=S_min+double(j1)/double(NNODE_1D-1)*S_range;
   s[2]=S_min+double(j2)/double(NNODE_1D-1)*S_range;
  }

 /// Get the local fraction of node j in the element
 void local_fraction_of_node(const unsigned &j,
                             Vector<double> &s_fraction)
  {
   s_fraction.resize(3);
   unsigned j0 = j%NNODE_1D;
   unsigned j1=(j/NNODE_1D)%NNODE_1D;
   unsigned j2=j/(NNODE_1D*NNODE_1D);
   s_fraction[0]= double(j0)/double(NNODE_1D-1);
   s_fraction[1]= double(j1)/double(NNODE_1D-1);
   s_fraction[2]= double(j2)/double(NNODE_1D-1);
  }

 /// This function returns the local fraction of any nodes in the n-th positoin
 /// in a one dimensional expansion along the i-th local coordinate
 inline double local_one_d_fraction_of_node(const unsigned &n1d,
                                            const unsigned &i)
  {
   //It's just the value of the node divided by the number of 1-D nodes
   return double(n1d)/double(NNODE_1D-1);
  }

 /// Get the node at the specified local coordinate
 Node* get_node_at_local_coordinate(const Vector<double> &s) const;
 
 /// Number of nodes along each element edge
 unsigned nnode_1d() const {return NNODE_1D;}

 /// \short Return the number of actual plot points for paraview
 /// plot with parameter nplot.
 unsigned nplot_points_paraview(const unsigned& nplot) const
  {
   return nplot*nplot*nplot;
  }

 /// \short Return the number of local sub-elements for paraview plot with 
 /// parameter nplot.
 unsigned nsub_elements_paraview(const unsigned& nplot) const 
  {
   return (nplot-1)*(nplot-1)*(nplot-1);
  }

 /// \short Fill in the offset information for paraview plot.
 /// Needs to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 void write_paraview_output_offset_information(std::ofstream& file_out,
                                               const unsigned& nplot,
                                               unsigned& counter) const
  {
   // Number of local elements we want to plot over
   unsigned plot=nsub_elements_paraview(nplot);
   
   for(unsigned j=0;j<plot;j+=(nplot-1)*(nplot-1)+1)
    {
     // To keep track of how many cross-sections we've looped over
     unsigned r=((j-(j%((nplot-1)*(nplot-1))))/((nplot-1)*(nplot-1)));

     // loop over all the elemnets in this sublevel
     unsigned sub_plot=(nplot-1)*(nplot-1);

     // loops over the i-th local element in parent element
     for(unsigned i=0;i<sub_plot;i++)
      {
       unsigned d=((i-(i%(nplot-1)))/(nplot-1));
       
     
       // Lower level of rectangle
       file_out << i%(nplot-1)+d*nplot+r*nplot*nplot+counter << " "
                << i%(nplot-1)+1+d*nplot+r*nplot*nplot+counter << " "
                << i%(nplot-1)+1+(d+1)*nplot+r*nplot*nplot+counter << " "
                << i%(nplot-1)+(d+1)*nplot+r*nplot*nplot+counter << " "
        
       // Upper level of rectangle
                << i%(nplot-1)+d*nplot+(r+1)*nplot*nplot+counter << " "
                << i%(nplot-1)+1+d*nplot+(r+1)*nplot*nplot+counter << " "
                << i%(nplot-1)+1+(d+1)*nplot+(r+1)*nplot*nplot+counter << " "
                << i%(nplot-1)+(d+1)*nplot+(r+1)*nplot*nplot+counter
                << std::endl;
      }
    }
   counter+=nplot_points_paraview(nplot);
  }

 /// \short Return the paraview element type.
 /// Needs to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 /// Use type "VTK_HEXAHEDRON" (== 12) for 2D quad elements
 void write_paraview_type(std::ofstream& file_out,
                          const unsigned& nplot) const
  {
   unsigned local_loop=nsub_elements_paraview(nplot);
   for(unsigned i=0;i<local_loop;i++)
    {
     file_out << "12" << std::endl;
    }
  }

 /// \short Return the offsets for the paraview sub-elements. Needs 
 /// to be implemented for each new geometric element type; see
 /// http://www.vtk.org/VTK/img/file-formats.pdf
 void write_paraview_offsets(std::ofstream& file_out,
                             const unsigned& nplot,
                             unsigned& offset_sum) const
  {
   unsigned local_loop=nsub_elements_paraview(nplot);
   for(unsigned i=0;i<local_loop;i++)
    {
     offset_sum+=8;
     file_out << offset_sum << std::endl;
    }
  }

 /// Output
 void output(std::ostream &outfile);

 /// Output at n_plot points
 void output(std::ostream &outfile, const unsigned &n_plot);

 /// C-style output 
 void output(FILE* file_pt);

 /// C_style output at n_plot points
 void output(FILE* file_pt, const unsigned &n_plot);

 /// \short  Get cector of local coordinates of plot point i (when plotting 
 /// nplot points in each "coordinate direction).
 void get_s_plot(
  const unsigned& i,
  const unsigned& nplot,
  Vector<double>& s,
  const bool& use_equally_spaced_interior_sample_points=false) const
 {
  if (nplot>1)
   {
    unsigned i01=i%(nplot*nplot);
    unsigned i0=i01%nplot;
    unsigned i1=(i01-i0)/nplot;
    unsigned i2=(i-i01)/(nplot*nplot);
    
    s[0]=-1.0+2.0*double(i0)/double(nplot-1);
    s[1]=-1.0+2.0*double(i1)/double(nplot-1);
    s[2]=-1.0+2.0*double(i2)/double(nplot-1); 
    if (use_equally_spaced_interior_sample_points)
     {
      double range=2.0;
      double dx_new=range/double(nplot);
      double range_new=double(nplot-1)*dx_new;
      s[0]=-1.0+0.5*dx_new+range_new*(1.0+s[0])/range;
      s[1]=-1.0+0.5*dx_new+range_new*(1.0+s[1])/range;
      s[2]=-1.0+0.5*dx_new+range_new*(1.0+s[2])/range;
     }
   }
  else
   {
    s[0]=0.0;
    s[1]=0.0;
    s[2]=0.0;
   }
 }

 /// \short Return string for tecplot zone header (when plotting 
 /// nplot points in each "coordinate direction)
 std::string tecplot_zone_string(const unsigned& nplot) const
  {
   std::ostringstream header;
   header << "ZONE I=" << nplot << ", J=" << nplot << ", K=" 
          << nplot << "\n";
   return header.str();
  }

 /// Return total number of plot points (when plotting 
 /// nplot points in each "coordinate direction)
 unsigned nplot_points(const unsigned& nplot) const
  {
   unsigned DIM=3;
   unsigned np=1;
   for (unsigned i=0;i<DIM;i++) {np*=nplot;}
   return np;
  };

 /// Get the number of the ith node on face face_index in the bulk node
 /// vector. 
 unsigned get_bulk_node_number(const int& face_index,
                               const unsigned& i) const
 {
  face_node_number_error_check(i);

  const unsigned nn1d = nnode_1d();

  if(face_index == -1) {return i*nn1d;}
  else if(face_index == +1) {return i*nn1d + (nn1d-1);}
  else if(face_index == -2) {return i%nn1d + (i/nn1d)*nn1d*nn1d;}
  else if(face_index == +2) {return i%nn1d + (i/nn1d)*nn1d*nn1d + (nn1d*(nn1d-1));} 
  else if(face_index == -3) {return i;} 
  else if(face_index == +3) {return i+(nn1d*nn1d)*(nn1d-1);}
  else
   {
    std::string err = "Face index should be in {-1, +1, -2, +2, -3, +3}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   }
 }

 /// Get the sign of the outer unit normal on the face given by face_index.
 int face_outer_unit_normal_sign(const int& face_index) const
 {
  if(face_index == -3) {return -1;} 
  else if(face_index == +3) {return 1;}
  else if(face_index == -2) {return 1;}
  else if(face_index == 2) {return -1;}
  else if(face_index == -1) {return -1;}
  else if(face_index == 1) {return 1;}
  else
   {
    std::string err = "Face index should be in {-1, +1, -2, +2, -3, +3}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   }
 }

 /// Get a pointer to the function mapping face coordinates to bulk coordinates
 CoordinateMappingFctPt face_to_bulk_coordinate_fct_pt
 (const int& face_index) const
 {
  if(face_index == 1) {return &QElement3FaceToBulkCoordinates::face3;}
  else if (face_index == -1) {return &QElement3FaceToBulkCoordinates::face0;}
  else if (face_index == -2) {return &QElement3FaceToBulkCoordinates::face1;}
  else if (face_index == 2) {return &QElement3FaceToBulkCoordinates::face4;}
  else if (face_index == -3) {return &QElement3FaceToBulkCoordinates::face2;}
  else if (face_index == 3) {return &QElement3FaceToBulkCoordinates::face5;}
  else
   {
    std::string err = "Face index should be in {-1, +1}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   } 
 }

 /// Get a pointer to the derivative of the mapping from face to bulk
 /// coordinates.
 BulkCoordinateDerivativesFctPt bulk_coordinate_derivatives_fct_pt
 (const int& face_index) const
 {
  if(face_index == 1) {return &QElement3BulkCoordinateDerivatives::faces0;}
  else if (face_index == -1) {return &QElement3BulkCoordinateDerivatives::faces0;}
  else if (face_index == -2) {return &QElement3BulkCoordinateDerivatives::faces1;}
  else if (face_index == 2) {return &QElement3BulkCoordinateDerivatives::faces1;}
  else if (face_index == -3) {return &QElement3BulkCoordinateDerivatives::faces2;}
  else if (face_index == 3) {return &QElement3BulkCoordinateDerivatives::faces2;}
  else
   {
    std::string err = "Face index should be in {-1, +1}.";
    throw OomphLibError(err, OOMPH_EXCEPTION_LOCATION,
                        OOMPH_CURRENT_FUNCTION);
   } 
 }

};
  
//=======================================================================
/// SolidQElement elements are quadrilateral elements whose 
/// derivatives also include those based upon the lagrangian 
/// positions of the nodes.
/// They are the basis for solid mechanics elements
//=======================================================================
template <unsigned DIM, unsigned NNODE_1D> 
class SolidQElement
{};


//=======================================================================
///SolidQElement elements, specialised to one spatial dimension
//=======================================================================
template <unsigned NNODE_1D>
class SolidQElement<1,NNODE_1D> : public virtual QElement<1,NNODE_1D>, 
 public virtual QSolidElementBase
{
  public:

 /// Constructor
 SolidQElement() : QElement<1,NNODE_1D>(), SolidFiniteElement() 
  {
   //Set the Lagrangian dimension of the element
   set_lagrangian_dimension(1);
  }
 
 /// Broken copy constructor
 SolidQElement(const SolidQElement&) 
  { 
   BrokenCopy::broken_copy("SolidQElement");
  } 
 
 /// Broken assignment operator
 /*void operator=(const SolidQElement&) 
  {
   BrokenCopy::broken_assign("SolidQElement");
   }*/

///Overload the output function
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}

 /// Output at n_plot points
 inline void output(std::ostream &outfile, const unsigned &n_plot);

 /// C-style output 
 void output(FILE* file_pt) {FiniteElement::output(file_pt);}

 /// C_style output at n_plot points
 inline void output(FILE* file_pt, const unsigned &n_plot);

 /// \short Build the lower-dimensional FaceElement (an element of type
 /// SolidQElement<0,NNODE_1D>).
 /// The face index takes one of two values corresponding
 /// to the two possible faces:
 /// -1 (Left)  s[0] = -1.0
 /// +1 (Right) s[0] =  1.0
 inline void build_face_element(const int &face_index, 
                                FaceElement* face_element_pt);


};

//For the dumb Intel 9.0 compiler, these need to live in here

///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////
// SolidQElements
///////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////
// 1D SolidQElements
///////////////////////////////////////////////////////////////////////////


//=======================================================================
/// The output function for n_plot points in each coordinate direction
/// for the 1D element
//=======================================================================
template <unsigned NNODE_1D>
void SolidQElement<1,NNODE_1D>::output(std::ostream &outfile, 
                                       const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(1);

 //Tecplot header info 
 outfile << "ZONE I=" << n_plot << std::endl;

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();
 
 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();

 //Loop over element nodes
 for(unsigned l=0;l<n_plot;l++)
  {
   s[0] = -1.0 + l*2.0/(n_plot-1);

   //Output the Eulerian coordinates
   for(unsigned i=0;i<n_dim;i++)
    {
     outfile << QElement<1,NNODE_1D>::interpolated_x(s,i) << " " ;
    }
   //Output the Lagranian coordinates
   for(unsigned i=0;i<n_lagr;i++)
    {
     outfile << SolidQElement<1,NNODE_1D>::interpolated_xi(s,i) << " " ;
    }
   outfile << std::endl;
  }
  outfile << std::endl;
}

//=======================================================================
/// C-style output function for n_plot points in each coordinate direction
/// for the 1D element
//=======================================================================
template <unsigned NNODE_1D>
void SolidQElement<1,NNODE_1D>::output(FILE* file_pt,
                                      const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(1);

 //Tecplot header info 
 fprintf(file_pt,"ZONE I=%i\n",n_plot);

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();
 
 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();

 //Loop over element nodes
 for(unsigned l=0;l<n_plot;l++)
  {
   s[0] = -1.0 + l*2.0/(n_plot-1);

   //Output the Eulerian coordinates
   for(unsigned i=0;i<n_dim;i++)
    {
     fprintf(file_pt,"%g ",QElement<1,NNODE_1D>::interpolated_x(s,i));
    }
   //Output the Lagranian coordinates
   for(unsigned i=0;i<n_lagr;i++)
    {
     fprintf(file_pt,"%g ",SolidQElement<1,NNODE_1D>::interpolated_xi(s,i));
    }
   fprintf(file_pt,"\n");
  }
 fprintf(file_pt,"\n");

}


//===========================================================
/// Function to setup geometrical information for lower-dimensional 
/// FaceElements (which are of type SolidQElement<0,1>).
//===========================================================
template<unsigned NNODE_1D>
void SolidQElement<1,NNODE_1D>::
build_face_element(const int &face_index,
                   FaceElement *face_element_pt)
{
 //Build the standard non-solid FaceElement
 QElement<1,NNODE_1D>::build_face_element(face_index,face_element_pt);

 //Set the Lagrangian dimension from the first node of the present element
 dynamic_cast<SolidFiniteElement*>(face_element_pt)->
  set_lagrangian_dimension(static_cast<SolidNode*>(node_pt(0))->nlagrangian());
}
 


//=======================================================================
/// SolidQElement elements, specialised to two spatial dimensions
//=======================================================================
template <unsigned NNODE_1D>
class SolidQElement<2,NNODE_1D> : public virtual QElement<2,NNODE_1D>, 
 public virtual QSolidElementBase
{
  public:

 /// Constructor
 SolidQElement() : QElementBase(), QElement<2,NNODE_1D>(),
  SolidFiniteElement(), QSolidElementBase()
  {
   //Set the Lagrangian dimension of the element
   set_lagrangian_dimension(2);
  }

 /// Broken copy constructor
 SolidQElement(const SolidQElement&) 
  { 
   BrokenCopy::broken_copy("SolidQElement");
  } 
 
 /// Broken assignment operator
 /*void operator=(const SolidQElement&) 
  {
   BrokenCopy::broken_assign("SolidQElement");
   }*/

 /// Overload the output function
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}

 /// Output at n_plot^2 points
 inline void output(std::ostream &outfile, const unsigned &n_plot);

 /// C-style output
 void output(FILE* file_pt){FiniteElement::output(file_pt);}

 /// C_style output at n_plot points
 inline void output(FILE* file_pt, const unsigned &n_plot);


 /// \short Build the lower-dimensional FaceElement (an element of type
 /// SolidQElement<1,NNODE_1D>).The face index takes one of four values
 /// corresponding to the four possible faces:
 /// -1 (West)  s[0] = -1.0
 /// +1 (East)  s[0] =  1.0
 /// -2 (South) s[1] = -1.0
 /// +2 (North) s[1] =  1.0
 inline void build_face_element(const int &face_index,
                                FaceElement* face_element_pt);

};





///////////////////////////////////////////////////////////////////////////
// 2D SolidQElements
///////////////////////////////////////////////////////////////////////////

//===========================================================
/// The output function for any number of points per element
//===========================================================
template <unsigned NNODE_1D>
void SolidQElement<2,NNODE_1D>::output(std::ostream &outfile, 
                                       const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(2);

 //Tecplot header info 
 outfile << "ZONE I=" << n_plot << ", J=" << n_plot << std::endl;

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();
 
 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();

 //Loop over plot points
 for(unsigned l2=0;l2<n_plot;l2++)
  {
   s[1] = -1.0 + l2*2.0/(n_plot-1);
   for(unsigned l1=0;l1<n_plot;l1++)
    {
     s[0] = -1.0 + l1*2.0/(n_plot-1);
     
     //Output the Eulerian coordinates
     for(unsigned i=0;i<n_dim;i++)
      {
       outfile << QElement<2,NNODE_1D>::interpolated_x(s,i) << " " ;
      }
     //Output the Lagranian coordinates
     for(unsigned i=0;i<n_lagr;i++)
      {
       outfile << SolidQElement<2,NNODE_1D>::interpolated_xi(s,i) << " " ;
      }
     outfile << std::endl;
    }
  }
 outfile << std::endl;
}




//====================================================================
/// C-style output function for any number of points per element
//====================================================================
template <unsigned NNODE_1D>
void SolidQElement<2,NNODE_1D>::output(FILE* file_pt,
                                       const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(2);

 //Tecplot header info 
 fprintf(file_pt,"ZONE I=%i, J=%i\n",n_plot,n_plot);

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();
 
 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();

 //Loop over plot points
 for(unsigned l2=0;l2<n_plot;l2++)
  {
   s[1] = -1.0 + l2*2.0/(n_plot-1);
   for(unsigned l1=0;l1<n_plot;l1++)
    {
     s[0] = -1.0 + l1*2.0/(n_plot-1);
     
     //Output the Eulerian coordinates
     for(unsigned i=0;i<n_dim;i++)
      {
       fprintf(file_pt,"%g ",QElement<2,NNODE_1D>::interpolated_x(s,i));
      }
     //Output the Lagranian coordinates
     for(unsigned i=0;i<n_lagr;i++)
      {
       fprintf(file_pt,"%g ",SolidQElement<2,NNODE_1D>::interpolated_xi(s,i));
      }
     fprintf(file_pt,"\n");
    }
  }
 fprintf(file_pt,"\n");
}


//===========================================================
/// Function to setup geometrical information for lower-dimensional 
/// FaceElements (which are of type SolidQElement<1,NNODE_1D>).
//===========================================================
template<unsigned NNODE_1D>
void SolidQElement<2,NNODE_1D>::
build_face_element(const int &face_index, FaceElement *face_element_pt)
{
 //Build the standard non-solid FaceElement
 QElement<2,NNODE_1D>::build_face_element(face_index,face_element_pt);
 
 //Set the Lagrangian dimension from the first node of the present element
 dynamic_cast<SolidFiniteElement*>(face_element_pt)->
  set_lagrangian_dimension(static_cast<SolidNode*>(node_pt(0))->nlagrangian());
}



//=======================================================================
/// SolidQElement elements, specialised to three spatial dimensions
//=======================================================================
template <unsigned NNODE_1D>
class SolidQElement<3,NNODE_1D> : public virtual QElement<3,NNODE_1D>, 
 public virtual QSolidElementBase
{

  public:

 /// Constructor
 SolidQElement() : QElementBase(), QElement<3,NNODE_1D>(),
  SolidFiniteElement(), QSolidElementBase()
  {
   //Set the Lagrangian dimension of the element
   set_lagrangian_dimension(3);
  }

 /// Broken copy constructor
 SolidQElement(const SolidQElement&) 
  { 
   BrokenCopy::broken_copy("SolidQElement");
  } 
 
 /// Broken assignment operator
 /*void operator=(const SolidQElement&) 
  {
   BrokenCopy::broken_assign("SolidQElement");
   }*/

 /// Overload the output function
 void output(std::ostream &outfile) {FiniteElement::output(outfile);}

 /// Output at n_plot^2 points
 inline void output(std::ostream &outfile, const unsigned &n_plot);

 /// C-style output
 void output(FILE* file_pt){FiniteElement::output(file_pt);}

 /// C_style output at n_plot points
 inline void output(FILE* file_pt, const unsigned &n_plot);


 /// \short Build the lower-dimensional FaceElement (an element of type
 /// SolidQElement<2,NNODE_1D>). The face index takes of one 
 /// six values corresponding
 /// to the six possible faces:
 /// -1 (Left)   s[0] = -1.0
 /// +1 (Right)  s[0] =  1.0
 /// -2 (Down)   s[1] = -1.0
 /// +2 (Up)     s[1] =  1.0
 /// -3 (Back)   s[2] = -1.0
 /// +3 (Front)  s[2] =  1.0
 inline void build_face_element(const int &face_index,
                                FaceElement* face_element_pt);

};





///////////////////////////////////////////////////////////////////////////
// 3D SolidQElements
///////////////////////////////////////////////////////////////////////////

//===========================================================
/// The output function for any number of points per element
//===========================================================
template <unsigned NNODE_1D>
void SolidQElement<3,NNODE_1D>::output(std::ostream &outfile, 
                                       const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(3);

 //Tecplot header info 
 outfile << "ZONE I=" << n_plot << ", J=" << n_plot 
  << ", K=" << n_plot << std::endl;

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();
 
 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();

 //Loop over plot points
 for(unsigned l3=0;l3<n_plot;l3++)
  {
   s[2] = -1.0 + l3*2.0/(n_plot-1);
   for(unsigned l2=0;l2<n_plot;l2++)
    {
     s[1] = -1.0 + l2*2.0/(n_plot-1);
     for(unsigned l1=0;l1<n_plot;l1++)
      {
       s[0] = -1.0 + l1*2.0/(n_plot-1);
       
       //Output the Eulerian coordinates
       for(unsigned i=0;i<n_dim;i++)
        {
         outfile << QElement<3,NNODE_1D>::interpolated_x(s,i) << " " ;
        }
       //Output the Lagranian coordinates
       for(unsigned i=0;i<n_lagr;i++)
        {
         outfile << SolidQElement<3,NNODE_1D>::interpolated_xi(s,i) << " " ;
        }
       outfile << std::endl;
      }
    }
  }
 outfile << std::endl;
}


//====================================================================
/// C-style output function for any number of points per element
//====================================================================
template <unsigned NNODE_1D>
void SolidQElement<3,NNODE_1D>::output(FILE* file_pt,
                                       const unsigned &n_plot)
{
 //Local variables
 Vector<double> s(3);

 //Tecplot header info 
 fprintf(file_pt,"ZONE I=%i, J=%i, K=%i\n",n_plot,n_plot,n_plot);

 //Find the dimension of the nodes
 unsigned n_dim = this->nodal_dimension();
 
 //Find the Lagrangian dimension of the first node
 unsigned n_lagr = static_cast<SolidNode*>(node_pt(0))->nlagrangian();

 //Loop over plot points
 for(unsigned l3=0;l3<n_plot;l3++)
  {
   s[2] = -1.0 + l3*2.0/(n_plot-1);
   for(unsigned l2=0;l2<n_plot;l2++)
    {
     s[1] = -1.0 + l2*2.0/(n_plot-1);
     for(unsigned l1=0;l1<n_plot;l1++)
      {
       s[0] = -1.0 + l1*2.0/(n_plot-1);
       
       //Output the Eulerian coordinates
       for(unsigned i=0;i<n_dim;i++)
        {
         fprintf(file_pt,"%g ",QElement<3,NNODE_1D>::interpolated_x(s,i));
        }
       //Output the Lagranian coordinates
       for(unsigned i=0;i<n_lagr;i++)
        {
         fprintf(file_pt,"%g ",
                 SolidQElement<3,NNODE_1D>::interpolated_xi(s,i));
        }
       fprintf(file_pt,"\n");
      }
    }
  }
 fprintf(file_pt,"\n");
}


//===========================================================
/// Function to setup geometrical information for lower-dimensional 
/// FaceElements (which are of type SolidQElement<1,NNODE_1D>).
//===========================================================
template<unsigned NNODE_1D>
void SolidQElement<3,NNODE_1D>::
 build_face_element(const int &face_index,
                    FaceElement *face_element_pt)
{
 //Build the standard non-solid FaceElement
 QElement<3,NNODE_1D>::build_face_element(face_index,face_element_pt);
 
 //Set the Lagrangian dimension from the first node of the present element
 dynamic_cast<SolidFiniteElement*>(face_element_pt)->
  set_lagrangian_dimension(static_cast<SolidNode*>(node_pt(0))->nlagrangian());
}


//==============================================================
/// A class that is used to template the refineable Q elements
/// by dimension. It's really nothing more than a policy class
//=============================================================
template<unsigned DIM>
class RefineableQElement
{
  public:

 /// Empty constuctor
 RefineableQElement() {}
};

//==============================================================
/// A class that is used to template the p-refineable Q elements
/// by dimension. It's really nothing more than a policy class.
/// The default template parameter ensures that these elements
/// inherit from the QElement of the correct type if they start
/// with a p-order higher than linear (e.g. Navier-Stokes Elements).
//=============================================================
template<unsigned DIM, unsigned INITIAL_NNODE_1D=2>
class PRefineableQElement
{
  public:

 /// Empty constuctor
 PRefineableQElement() {}
};

//==============================================================
/// A class that is used to template the solid refineable Q elements
/// by dimension. It's really nothing more than a policy class
//=============================================================
template<unsigned DIM>
class RefineableSolidQElement
{
  public:

 /// Empty constuctor
 RefineableSolidQElement() {}
};


}

#endif