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-2022 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
{
 
  //========================================================================
  //========================================================================
  class QElementGeometricBase : public virtual FiniteElement
  {
  public:
    QElementGeometricBase() {}
 
    QElementGeometricBase(const QElementGeometricBase&) = delete;
 
    // 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;*/
  };
 
 
 
 
  //========================================================================
  //========================================================================
  class QElementBase : public virtual QElementGeometricBase
  {
  public:
    QElementBase() : S_macro_ll_pt(0), S_macro_ur_pt(0) {}
 
    QElementBase(const QElementBase&) = delete;
 
    /*void operator=(const QElementBase&) = delete;*/
 
    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;
    }
 
    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;
    }
 
    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();
      }
    }
 
 
    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;
      }
 
      FiniteElement::set_macro_elem_pt(macro_elem_pt);
    }
 
 
    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];
    }
 
 
    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];
    }
 
 
    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];
    }
 
 
    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];
    }
 
    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);
    }
 
    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);
    }
 
    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)));
    }
 
    ElementGeometry::ElementGeometry element_geometry() const
    {
      return ElementGeometry::Q;
    }
 
  private:
    Vector<double>* S_macro_ll_pt;
 
    Vector<double>* S_macro_ur_pt;
  };
 
 
 
 
  //========================================================================
  //========================================================================
  class QSolidElementBase : public virtual QElementBase,
                            public virtual SolidFiniteElement
  {
  public:
    QSolidElementBase(){};
 
    QSolidElementBase(const QSolidElementBase&) = delete;
 
    /*void operator=(const QSolidElementBase&) = delete;*/
 
    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);
    }
 
    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);
    }
 
    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);
      }
    }
  };
 
 
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class QElement
  {
  };
 
  //=======================================================================
  //=======================================================================
  class LineElementBase : public virtual QElementBase
  {
  public:
    LineElementBase() {}
 
    virtual unsigned nvertex_node() const = 0;
 
    virtual Node* vertex_node_pt(const unsigned& j) const = 0;
  };
 
  //=======================================================================
  //=======================================================================
  template<unsigned NNODE_1D>
  class QElement<1, NNODE_1D> : public virtual LineElementBase
  {
  private:
    // 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:
    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);
    }
 
    QElement(const QElement&) = delete;
 
    /*void operator=(const QElement&) = delete;*/
 
    void shape(const Vector<double>& s, Shape& psi) const;
 
    void dshape_local(const Vector<double>& s,
                      Shape& psi,
                      DShape& dpsids) const;
 
    void d2shape_local(const Vector<double>& s,
                       Shape& psi,
                       DShape& dpsids,
                       DShape& d2psids) const;
 
    double invert_jacobian_mapping(const DenseMatrix<double>& jacobian,
                                   DenseMatrix<double>& inverse_jacobian) const
    {
      return FiniteElement::invert_jacobian<1>(jacobian, inverse_jacobian);
    }
 
    double s_min() const
    {
      return -1.0;
    }
 
    double s_max() const
    {
      return 1.0;
    }
 
    unsigned nvertex_node() const
    {
      return 2;
    }
 
    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;
    }
 
    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());
    }
 
 
    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);
    }
 
 
    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);
    }
 
    Node* get_node_at_local_coordinate(const Vector<double>& s) const;
 
    unsigned nnode_1d() const
    {
      return NNODE_1D;
    }
 
    unsigned nplot_points_paraview(const unsigned& nplot) const
    {
      return nplot;
    }
 
    unsigned nsub_elements_paraview(const unsigned& nplot) const
    {
      return (nplot - 1);
    }
 
    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);
    }
 
    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;
      }
    }
 
    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;
      }
    }
 
    void output(std::ostream& outfile);
 
    void output(std::ostream& outfile, const unsigned& n_plot);
 
    void output(FILE* file_pt);
 
    void output(FILE* file_pt, const unsigned& n_plot);
 
 
    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;
      }
    }
 
    std::string tecplot_zone_string(const unsigned& nplot) const
    {
      std::ostringstream header;
      header << "ZONE I=" << nplot << "\n";
      return header.str();
    }
 
    unsigned nplot_points(const unsigned& nplot) const
    {
      unsigned DIM = 1;
      unsigned np = 1;
      for (unsigned i = 0; i < DIM; i++)
      {
        np *= nplot;
      }
      return np;
    }
 
    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);
      }
    }
 
    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;
    }
 
    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);
      }
    }
 
    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);
      }
    }
  };
 
  //=======================================================================
  //=======================================================================
  class QuadElementBase : public virtual QElementBase
  {
  public:
    QuadElementBase() {}
 
    virtual unsigned nvertex_node() const = 0;
 
    virtual Node* vertex_node_pt(const unsigned& j) const = 0;
  };
 
  //=======================================================================
  //=======================================================================
  template<unsigned NNODE_1D>
  class QElement<2, NNODE_1D> : public virtual QuadElementBase
  {
  private:
    // 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:
    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);
    }
 
    QElement(const QElement&) = delete;
 
    /*void operator=(const QElement&) = delete;*/
 
    void shape(const Vector<double>& s, Shape& psi) const;
 
    void dshape_local(const Vector<double>& s,
                      Shape& psi,
                      DShape& dpsids) const;
 
    void d2shape_local(const Vector<double>& s,
                       Shape& psi,
                       DShape& dpsids,
                       DShape& d2psids) const;
 
    double invert_jacobian_mapping(const DenseMatrix<double>& jacobian,
                                   DenseMatrix<double>& inverse_jacobian) const
    {
      return FiniteElement::invert_jacobian<2>(jacobian, inverse_jacobian);
    }
 
    double s_min() const
    {
      return -1.0;
    }
 
    double s_max() const
    {
      return 1.0;
    }
 
 
    unsigned nvertex_node() const
    {
      return 4;
    }
 
    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;
    }
 
 
    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;
    }
 
    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);
    }
 
    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);
    }
 
    Node* get_node_at_local_coordinate(const Vector<double>& s) const;
 
    unsigned nnode_1d() const
    {
      return NNODE_1D;
    }
 
    unsigned nplot_points_paraview(const unsigned& nplot) const
    {
      return nplot * nplot;
    }
 
    unsigned nsub_elements_paraview(const unsigned& nplot) const
    {
      return (nplot - 1) * (nplot - 1);
    }
 
    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);
    }
 
    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;
      }
    }
 
    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;
      }
    }
 
    void output(std::ostream& outfile);
 
    void output(std::ostream& outfile, const unsigned& n_plot);
 
    void output(FILE* file_pt);
 
    void output(FILE* file_pt, const unsigned& n_plot);
 
 
    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;
      }
    }
 
    std::string tecplot_zone_string(const unsigned& nplot) const
    {
      std::ostringstream header;
      header << "ZONE I=" << nplot << ", J=" << nplot << "\n";
      return header.str();
    }
 
    unsigned nplot_points(const unsigned& nplot) const
    {
      unsigned DIM = 2;
      unsigned np = 1;
      for (unsigned i = 0; i < DIM; i++)
      {
        np *= nplot;
      }
      return np;
    };
 
    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);
      }
    }
 
    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);
      }
    }
 
    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);
      }
    }
 
    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);
      }
    }
  };
 
  //=======================================================================
  //=======================================================================
  class BrickElementBase : public virtual QElementBase
  {
  public:
    BrickElementBase() {}
 
    virtual unsigned nvertex_node() const = 0;
 
    virtual Node* vertex_node_pt(const unsigned& j) const = 0;
  };
 
  //=======================================================================
  //=======================================================================
  template<unsigned NNODE_1D>
  class QElement<3, NNODE_1D> : public virtual BrickElementBase
  {
  private:
    // 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:
    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);
    }
 
 
    QElement(const QElement&) = delete;
 
    /*void operator=(const QElement&) = delete;*/
 
    void shape(const Vector<double>& s, Shape& psi) const;
 
    void dshape_local(const Vector<double>& s,
                      Shape& psi,
                      DShape& dpsids) const;
 
    void d2shape_local(const Vector<double>& s,
                       Shape& psi,
                       DShape& dpsids,
                       DShape& d2psids) const;
 
 
    double invert_jacobian_mapping(const DenseMatrix<double>& jacobian,
                                   DenseMatrix<double>& inverse_jacobian) const
    {
      return FiniteElement::invert_jacobian<3>(jacobian, inverse_jacobian);
    }
 
    double s_min() const
    {
      return -1.0;
    }
 
    double s_max() const
    {
      return 1.0;
    }
 
    unsigned nvertex_node() const
    {
      return 8;
    }
 
    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;
    }
 
    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;
    }
 
    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);
    }
 
    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);
    }
 
    Node* get_node_at_local_coordinate(const Vector<double>& s) const;
 
    unsigned nnode_1d() const
    {
      return NNODE_1D;
    }
 
    unsigned nplot_points_paraview(const unsigned& nplot) const
    {
      return nplot * nplot * nplot;
    }
 
    unsigned nsub_elements_paraview(const unsigned& nplot) const
    {
      return (nplot - 1) * (nplot - 1) * (nplot - 1);
    }
 
    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);
    }
 
    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;
      }
    }
 
    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;
      }
    }
 
    void output(std::ostream& outfile);
 
    void output(std::ostream& outfile, const unsigned& n_plot);
 
    void output(FILE* file_pt);
 
    void output(FILE* file_pt, const unsigned& n_plot);
 
    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;
      }
    }
 
    std::string tecplot_zone_string(const unsigned& nplot) const
    {
      std::ostringstream header;
      header << "ZONE I=" << nplot << ", J=" << nplot << ", K=" << nplot
             << "\n";
      return header.str();
    }
 
    unsigned nplot_points(const unsigned& nplot) const
    {
      unsigned DIM = 3;
      unsigned np = 1;
      for (unsigned i = 0; i < DIM; i++)
      {
        np *= nplot;
      }
      return np;
    };
 
    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);
      }
    }
 
    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);
      }
    }
 
    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);
      }
    }
 
    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);
      }
    }
  };
 
  //=======================================================================
  //=======================================================================
  template<unsigned DIM, unsigned NNODE_1D>
  class SolidQElement
  {
  };
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned NNODE_1D>
  class SolidQElement<1, NNODE_1D> : public virtual QElement<1, NNODE_1D>,
                                     public virtual QSolidElementBase
  {
  public:
    SolidQElement() : QElement<1, NNODE_1D>(), SolidFiniteElement()
    {
      // Set the Lagrangian dimension of the element
      set_lagrangian_dimension(1);
    }
 
    SolidQElement(const SolidQElement&) = delete;
 
    /*void operator=(const SolidQElement&) = delete;*/
 
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    inline void output(std::ostream& outfile, const unsigned& n_plot);
 
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    inline void output(FILE* file_pt, const unsigned& n_plot);
 
    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
 
 
  //=======================================================================
  //=======================================================================
  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;
  }
 
  //=======================================================================
  //=======================================================================
  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");
  }
 
 
  //===========================================================
  //===========================================================
  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());
  }
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned NNODE_1D>
  class SolidQElement<2, NNODE_1D> : public virtual QElement<2, NNODE_1D>,
                                     public virtual QSolidElementBase
  {
  public:
    SolidQElement()
      : QElementBase(),
        QElement<2, NNODE_1D>(),
        SolidFiniteElement(),
        QSolidElementBase()
    {
      // Set the Lagrangian dimension of the element
      set_lagrangian_dimension(2);
    }
 
    SolidQElement(const SolidQElement&) = delete;
 
    /*void operator=(const SolidQElement&) = delete;*/
 
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    inline void output(std::ostream& outfile, const unsigned& n_plot);
 
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    inline void output(FILE* file_pt, const unsigned& n_plot);
 
 
    inline void build_face_element(const int& face_index,
                                   FaceElement* face_element_pt);
  };
 
 
  // 2D SolidQElements
 
  //===========================================================
  //===========================================================
  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;
  }
 
 
  //====================================================================
  //====================================================================
  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");
  }
 
 
  //===========================================================
  //===========================================================
  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());
  }
 
 
  //=======================================================================
  //=======================================================================
  template<unsigned NNODE_1D>
  class SolidQElement<3, NNODE_1D> : public virtual QElement<3, NNODE_1D>,
                                     public virtual QSolidElementBase
  {
  public:
    SolidQElement()
      : QElementBase(),
        QElement<3, NNODE_1D>(),
        SolidFiniteElement(),
        QSolidElementBase()
    {
      // Set the Lagrangian dimension of the element
      set_lagrangian_dimension(3);
    }
 
    SolidQElement(const SolidQElement&) = delete;
 
    /*void operator=(const SolidQElement&) = delete;*/
 
    void output(std::ostream& outfile)
    {
      FiniteElement::output(outfile);
    }
 
    inline void output(std::ostream& outfile, const unsigned& n_plot);
 
    void output(FILE* file_pt)
    {
      FiniteElement::output(file_pt);
    }
 
    inline void output(FILE* file_pt, const unsigned& n_plot);
 
 
    inline void build_face_element(const int& face_index,
                                   FaceElement* face_element_pt);
  };
 
 
  // 3D SolidQElements
 
  //===========================================================
  //===========================================================
  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;
  }
 
 
  //====================================================================
  //====================================================================
  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");
  }
 
 
  //===========================================================
  //===========================================================
  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());
  }
 
 
  //==============================================================
  //=============================================================
  template<unsigned DIM>
  class RefineableQElement
  {
  public:
    RefineableQElement() {}
  };
 
  //==============================================================
  //=============================================================
  template<unsigned DIM, unsigned INITIAL_NNODE_1D = 2>
  class PRefineableQElement
  {
  public:
    PRefineableQElement() {}
  };
 
  //==============================================================
  //=============================================================
  template<unsigned DIM>
  class RefineableSolidQElement
  {
  public:
    RefineableSolidQElement() {}
  };
 
 
} // namespace oomph
 
#endif