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// Copyright (C) 2006-2024 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&) = delete;
 
    /// 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&) = delete;*/
  };
 
 
  /// ///////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////
  /// ///////////////////////////////////////////////////////////////////
 
 
  //========================================================================
  /// 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&) = delete;
 
    /// Broken assignment operator
    /*void operator=(const QElementBase&) = delete;*/
 
    /// 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;
    }
 
    /// 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();
      }
    }
 
 
    /// 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);
    }
 
 
    /// 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];
    }
 
 
    /// 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];
    }
 
 
    /// 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];
    }
 
 
    /// 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];
    }
 
    /// 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);
    }
 
    /// 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&) = delete;
 
    /// Broken assignment operator
    /*void operator=(const QSolidElementBase&) = delete;*/
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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() {}
 
    /// Number of vertex nodes in the element
    virtual unsigned nvertex_node() const = 0;
 
    /// 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:
    /// 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&) = delete;
 
    /// Broken assignment operator
    /*void operator=(const QElement&) = delete;*/
 
    /// Calculate the geometric shape functions at local coordinate s
    void shape(const Vector<double>& s, Shape& psi) const;
 
    /// 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;
 
    /// 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;
 
    /// 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;
    }
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return 2;
    }
 
    /// 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;
    }
 
    /// Return the number of actual plot points for paraview
    /// plot with parameter nplot.
    unsigned nplot_points_paraview(const unsigned& nplot) const
    {
      return nplot;
    }
 
    /// Return the number of local sub-elements for paraview plot with
    /// parameter nplot.
    unsigned nsub_elements_paraview(const unsigned& nplot) const
    {
      return (nplot - 1);
    }
 
    /// 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);
    }
 
    /// 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;
      }
    }
 
    /// 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);
 
 
    ///  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;
      }
    }
 
    /// 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();
    }
 
    /// 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() {}
 
    /// Number of vertex nodes in the element
    virtual unsigned nvertex_node() const = 0;
 
    /// 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:
    /// 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&) = delete;
 
    /// Broken assignment operator
    /*void operator=(const QElement&) = delete;*/
 
    /// Calculate the geometric shape functions at local coordinate s
    void shape(const Vector<double>& s, Shape& psi) const;
 
    /// 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;
 
    /// 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;
 
    /// 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;
    }
 
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return 4;
    }
 
    /// 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;
    }
 
    /// Return the number of actual plot points for paraview
    /// plot with parameter nplot.
    unsigned nplot_points_paraview(const unsigned& nplot) const
    {
      return nplot * nplot;
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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;
      }
    }
 
    /// 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);
 
 
    ///  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;
      }
    }
 
    /// 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() {}
 
    /// Number of vertex nodes in the element
    virtual unsigned nvertex_node() const = 0;
 
    /// 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:
    /// 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&) = delete;
 
    /// Broken assignment operator
    /*void operator=(const QElement&) = delete;*/
 
    /// Calculate the geometric shape functions at local coordinate s
    void shape(const Vector<double>& s, Shape& psi) const;
 
    /// 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;
 
    /// 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;
 
 
    /// 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;
    }
 
    /// Number of vertex nodes in the element
    unsigned nvertex_node() const
    {
      return 8;
    }
 
    /// 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;
    }
 
    /// 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;
    }
 
    /// 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);
    }
 
    /// 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);
    }
 
    /// 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;
      }
    }
 
    /// 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);
 
    ///  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;
      }
    }
 
    /// 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&) = delete;
 
    /// Broken assignment operator
    /*void operator=(const SolidQElement&) = delete;*/
 
    /// 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);
 
    /// 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&) = delete;
 
    /// Broken assignment operator
    /*void operator=(const SolidQElement&) = delete;*/
 
    /// 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);
 
 
    /// 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&) = delete;
 
    /// Broken assignment operator
    /*void operator=(const SolidQElement&) = delete;*/
 
    /// 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);
 
 
    /// 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() {}
  };
 
 
} // namespace oomph
 
#endif