shape.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-2023 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//====================================================================
// This header file includes generic shape function classes
 
// Include guards to prevent multiple inclusions of the file
#ifndef OOMPH_SHAPE_HEADER
#define OOMPH_SHAPE_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 includes
#include "Vector.h"
#include "matrices.h"
#include "orthpoly.h"
 
namespace oomph
{
  //========================================================================
  /// A Class for shape functions. In simple cases, the shape functions
  /// have only one index that can be thought of as corresponding to the
  /// nodal points. In general, however, when quantities and
  /// their gradients are interpolated separately, the shape function have
  /// two indices: one corresponding to the nodal points, and the other
  /// to the "type" of quantity being interpolated: function, derivative, &c
  /// The second index can also represent the vector coordinate for
  /// vector-valued (Nedelec) shape functions.
  ///
  /// The implementation of Shape functions is designed to permit fast
  /// copying of entire sets of values by resetting the internal pointer
  /// to the data, Psi;
  /// functionality that is required, for example,
  ///  when setting the test functions
  /// in Galerkin elements and when reading pre-computed values of the shape
  /// functions.
  /// In general, we cannot know at construction time whether the pointer to
  /// the values will be reset or not and, therefore,
  /// whether the storage for values should be allocated by the object.
  /// We choose to allocate storage on construction and store an
  /// additional pointer Allocated_data that \b always addresses the storage
  /// allocated by the object. If the Psi pointer is reset then this storage
  /// will be "wasted", but only for the lifetime of the object. The cost for
  /// non-copied Shape functions is one additional pointer.
  //=========================================================================
  class Shape
  {
  protected:
    /// Pointer that addresses the storage that will be used to read and
    /// set the shape functions. The shape functions are packed into
    /// a flat array of doubles.
    double* Psi;
 
    /// Pointer that addresses the storage allocated by the object on
    /// construction. This will be the same as Psi if the object is not
    /// copied.
    double* Allocated_storage;
 
    /// Size of the first index of the shape function
    unsigned Index1;
 
    /// Size of the second index of the shape function
    unsigned Index2;
 
    /// Private function that checks whether the index is in range
    void range_check(const unsigned& i, const unsigned& j) const
    {
      // If an index is out of range, throw an error
      if ((i >= Index1) || (j >= Index2))
      {
        std::ostringstream error_stream;
        error_stream << "Range Error: ";
        if (i >= Index1)
        {
          error_stream << i << " is not in the range (0," << Index1 - 1 << ")"
                       << std::endl;
        }
        if (j >= Index2)
        {
          error_stream << j << " is not in the range (0," << Index2 - 1 << ")"
                       << std::endl;
        }
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
  public:
    /// Constructor for a single-index set of shape functions.
    Shape(const unsigned& N) : Index1(N), Index2(1)
    {
      Allocated_storage = new double[N];
      Psi = Allocated_storage;
    }
 
    /// Constructor for a two-index set of shape functions.
    Shape(const unsigned& N, const unsigned& M) : Index1(N), Index2(M)
    {
      Allocated_storage = new double[N * M];
      Psi = Allocated_storage;
    }
 
    /// Broken copy constructor
    Shape(const Shape& shape) = delete;
 
    /// Default constructor - just assigns a null pointers and zero index
    /// sizes.
    Shape() : Psi(0), Allocated_storage(0), Index1(0), Index2(0) {}
 
    /// The assignment operator does a shallow copy
    /// (resets the pointer to the data)
    void operator=(const Shape& shape)
    {
#ifdef PARANOID
      // Check the dimensions
      if ((shape.Index1 != Index1) || (shape.Index2 != Index2))
      {
        std::ostringstream error_stream;
        error_stream << "Cannot assign Shape object:" << std::endl
                     << "Indices do not match "
                     << "LHS: " << Index1 << " " << Index2
                     << ", RHS: " << shape.Index1 << " " << shape.Index2
                     << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      Psi = shape.Psi;
    }
 
    /// The assignment operator does a shallow copy
    /// (resets the pointer to the data)
    void operator=(Shape* const& shape_pt)
    {
#ifdef PARANOID
      // Check the dimensions
      if ((shape_pt->Index1 != Index1) || (shape_pt->Index2 != Index2))
      {
        std::ostringstream error_stream;
        error_stream << "Cannot assign Shape object:" << std::endl
                     << "Indices do not match "
                     << "LHS: " << Index1 << " " << Index2
                     << ", RHS: " << shape_pt->Index1 << " " << shape_pt->Index2
                     << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      Psi = shape_pt->Psi;
    }
 
    /// Destructor, clear up the memory allocated by the object
    ~Shape()
    {
      delete[] Allocated_storage;
      Allocated_storage = 0;
    }
 
    /// Change the size of the storage
    void resize(const unsigned& N, const unsigned& M = 1)
    {
      // Clear old storage
      delete[] Allocated_storage;
      Allocated_storage = 0;
      Psi = 0;
 
      // Allocate new storage
      Index1 = N;
      Index2 = M;
      Allocated_storage = new double[N * M];
      Psi = Allocated_storage;
    }
 
    /// Overload the bracket operator to provide access to values.
    inline double& operator[](const unsigned& i)
    {
#ifdef RANGE_CHECKING
      range_check(i, 0);
#endif
      return Psi[i * Index2];
    }
 
    /// Overload the bracket operator (const version)
    inline const double& operator[](const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      range_check(i, 0);
#endif
      return Psi[i * Index2];
    }
 
    /// Overload the round bracket operator to provide access to values.
    inline double& operator()(const unsigned& i)
    {
#ifdef RANGE_CHECKING
      range_check(i, 0);
#endif
      return Psi[i * Index2];
    }
 
    /// Overload the round bracket operator (const version)
    inline const double& operator()(const unsigned& i) const
    {
#ifdef RANGE_CHECKING
      range_check(i, 0);
#endif
      return Psi[i * Index2];
    }
 
    /// Overload the round bracket operator, allowing for two indices
    inline double& operator()(const unsigned& i, const unsigned& j)
    {
#ifdef RANGE_CHECKING
      range_check(i, j);
#endif
      return Psi[i * Index2 + j];
    }
 
    /// Overload the round bracket operator, allowing for two indices
    /// (const version)
    inline const double& operator()(const unsigned& i, const unsigned& j) const
    {
#ifdef RANGE_CHECKING
      range_check(i, j);
#endif
      return Psi[i * Index2 + j];
    }
 
    /// Return the range of index 1 of the shape function object
    inline unsigned nindex1() const
    {
      return Index1;
    }
 
    /// Return the range of index 2 of the shape function object
    inline unsigned nindex2() const
    {
      return Index2;
    }
  };
 
  //================================================================
  /// A Class for the derivatives of shape functions
  /// The class design is essentially the same as Shape, but there is
  /// on additional index that is used to indicate the coordinate direction in
  /// which the derivative is taken.
  //================================================================
  class DShape
  {
  private:
    /// Pointer that addresses the storage that will be used to read and
    /// set the shape-function derivatives. The values are packed into
    /// a flat array of doubles.
    double* DPsi;
 
    /// Pointer that addresses the storage allocated by the object on
    /// construction. This will be the same as DPsi if the object is not
    /// copied.
    double* Allocated_storage;
 
    /// Size of the first index of the shape function
    unsigned Index1;
 
    /// Size of the second index of the shape function
    unsigned Index2;
 
    /// Size of the third index of the shape function
    unsigned Index3;
 
    /// Private function that checks whether the indices are in range
    void range_check(const unsigned& i,
                     const unsigned& j,
                     const unsigned& k) const
    {
      // Check the first index
      if ((i >= Index1) || (j >= Index2) || (k >= Index3))
      {
        std::ostringstream error_stream;
        error_stream << "Range Error: ";
        if (i >= Index1)
        {
          error_stream << i << " is not in the range (0," << Index1 - 1 << ")"
                       << std::endl;
        }
        if (j >= Index2)
        {
          error_stream << j << " is not in the range (0," << Index2 - 1 << ")"
                       << std::endl;
        }
        if (k >= Index3)
        {
          error_stream << k << " is not in the range (0," << Index3 - 1 << ")"
                       << std::endl;
        }
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
    }
 
 
  public:
    /// Constructor with two parameters: a single-index shape function
    DShape(const unsigned& N, const unsigned& P)
      : Index1(N), Index2(1), Index3(P)
    {
      Allocated_storage = new double[N * P];
      DPsi = Allocated_storage;
    }
 
    /// Constructor with three paramters: a two-index shape function
    DShape(const unsigned& N, const unsigned& M, const unsigned& P)
      : Index1(N), Index2(M), Index3(P)
    {
      Allocated_storage = new double[N * M * P];
      DPsi = Allocated_storage;
    }
 
    /// Default constructor - just assigns a null pointers and zero index
    /// sizes.
    DShape() : DPsi(0), Allocated_storage(0), Index1(0), Index2(0), Index3(0) {}
 
    /// Broken copy constructor
    DShape(const DShape& dshape) = delete;
 
    /// The assignment operator does a shallow copy
    /// (resets the pointer to the data)
    void operator=(const DShape& dshape)
    {
#ifdef PARANOID
      // Check the dimensions
      if ((dshape.Index1 != Index1) || (dshape.Index2 != Index2) ||
          (dshape.Index3 != Index3))
      {
        std::ostringstream error_stream;
        error_stream << "Cannot assign DShape object:" << std::endl
                     << "Indices do not match "
                     << "LHS: " << Index1 << " " << Index2 << " " << Index3
                     << ", RHS: " << dshape.Index1 << " " << dshape.Index2
                     << " " << dshape.Index3 << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      DPsi = dshape.DPsi;
    }
 
    /// The assignment operator does a shallow copy
    /// (resets the pointer to the data)
    void operator=(DShape* const& dshape_pt)
    {
#ifdef PARANOID
      // Check the dimensions
      if ((dshape_pt->Index1 != Index1) || (dshape_pt->Index2 != Index2) ||
          (dshape_pt->Index3 != Index3))
      {
        std::ostringstream error_stream;
        error_stream << "Cannot assign Shape object:" << std::endl
                     << "Indices do not match "
                     << "LHS: " << Index1 << " " << Index2 << " " << Index3
                     << ", RHS: " << dshape_pt->Index1 << " "
                     << dshape_pt->Index2 << " " << dshape_pt->Index3
                     << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      DPsi = dshape_pt->DPsi;
    }
 
 
    /// Destructor, clean up the memory allocated by this object
    ~DShape()
    {
      delete[] Allocated_storage;
      Allocated_storage = 0;
    }
 
    /// Change the size of the storage. Note that (for some strange reason)
    /// index2 is the "optional" index, to conform with the existing
    /// constructor.
    void resize(const unsigned& N, const unsigned& P, const unsigned& M = 1)
    {
      // Clear old storage
      delete[] Allocated_storage;
      Allocated_storage = 0;
      DPsi = 0;
 
      // Allocate new storage
      Index1 = N;
      Index2 = M;
      Index3 = P;
      Allocated_storage = new double[N * M * P];
      DPsi = Allocated_storage;
    }
 
    /// Overload the round bracket operator for access to the data
    inline double& operator()(const unsigned& i, const unsigned& k)
    {
#ifdef RANGE_CHECKING
      range_check(i, 0, k);
#endif
      return DPsi[i * Index2 * Index3 + k];
    }
 
    /// Overload the round bracket operator (const version)
    inline const double& operator()(const unsigned& i, const unsigned& k) const
    {
#ifdef RANGE_CHECKING
      range_check(i, 0, k);
#endif
      return DPsi[i * Index2 * Index3 + k];
    }
 
    /// Overload the round bracket operator, with 3 indices
    inline double& operator()(const unsigned& i,
                              const unsigned& j,
                              const unsigned& k)
    {
#ifdef RANGE_CHECKING
      range_check(i, j, k);
#endif
      return DPsi[(i * Index2 + j) * Index3 + k];
    }
 
    /// Overload the round bracket operator (const version)
    inline const double& operator()(const unsigned& i,
                                    const unsigned& j,
                                    const unsigned& k) const
    {
#ifdef RANGE_CHECKING
      range_check(i, j, k);
#endif
      return DPsi[(i * Index2 + j) * Index3 + k];
    }
 
    /// Direct access to internal storage of data in flat-packed C-style
    /// column-major format. WARNING: Only for experienced users. Only
    /// use this if raw speed is of the essence, as in the solid mechanics
    /// problems.
    inline double& raw_direct_access(const unsigned long& i)
    {
      return DPsi[i];
    }
 
    /// Direct access to internal storage of data in flat-packed C-style
    /// column-major format. WARNING: Only for experienced users. Only
    /// use this if raw speed is of the essence, as in the solid mechanics
    /// problems.
    inline const double& raw_direct_access(const unsigned long& i) const
    {
      return DPsi[i];
    }
 
    /// Caculate the offset in flat-packed C-style, column-major format,
    /// required for a given i,j. WARNING: Only for experienced users. Only
    /// use this if raw speed is of the essence, as in the solid mechanics
    /// problems.
    unsigned offset(const unsigned long& i, const unsigned long& j) const
    {
      return (i * Index2 + j) * Index3 + 0;
    }
 
 
    /// Return the range of index 1 of the derivatives of the shape functions
    inline unsigned long nindex1() const
    {
      return Index1;
    }
 
    /// Return the range of index 2 of the derivatives of the shape functions
    inline unsigned long nindex2() const
    {
      return Index2;
    }
 
    /// Return the range of index 3 of the derivatives of the shape functions
    inline unsigned long nindex3() const
    {
      return Index3;
    }
  };
 
  //======================================================================
  /// A shape function with a deep copy constructor. This allows for use with
  /// stl
  ///  operations (e.g. manipulating vectors of shape functions). A seperate
  ///  class is needed because the basic shape function uses a shallow copy.
  //======================================================================
  class ShapeWithDeepCopy : public Shape
  {
  public:
    /// Constructor for a single-index set of shape functions.
    ShapeWithDeepCopy(const unsigned& N) : Shape(N) {}
 
    /// Constructor for a two-index set of shape functions.
    ShapeWithDeepCopy(const unsigned& N, const unsigned& M) : Shape(N, M) {}
 
    /// Default constructor
    ShapeWithDeepCopy() : Shape() {}
 
    /// Deep copy constructor
    ShapeWithDeepCopy(const ShapeWithDeepCopy& old_shape)
      : Shape(old_shape.Index1, old_shape.Index2)
    {
      for (unsigned i = 0; i < Index1 * Index2; i++)
      {
        Psi[i] = old_shape.Psi[i];
      }
    }
 
    /// Broken assignment operator
    void operator=(const ShapeWithDeepCopy& old_shape) = delete;
 
    /// Destructor, clear up the memory allocated by the object
    ~ShapeWithDeepCopy()
    {
      delete[] Allocated_storage;
      Allocated_storage = 0;
    }
  };
 
  /// /////////////////////////////////////////////////////////////////
  //
  // One dimensional shape functions and derivatives.
  // empty -- simply establishes the template parameters.
  //
  /// /////////////////////////////////////////////////////////////////
 
  namespace OneDimLagrange
  {
    /// Definition for 1D Lagrange shape functions. The
    /// value of all the shape functions at the local coordinate s
    /// are returned in the array Psi.
    template<unsigned NNODE_1D>
    void shape(const double& s, double* Psi)
    {
      std::ostringstream error_stream;
      error_stream << "One dimensional Lagrange shape functions "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Definition for derivatives of 1D Lagrange shape functions. The
    /// value of all the shape function derivatives at the local coordinate s
    /// are returned in the array DPsi.
    template<unsigned NNODE_1D>
    void dshape(const double& s, double* DPsi)
    {
      std::ostringstream error_stream;
      error_stream << "One dimensional Lagrange shape function derivatives "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Definition for second derivatives of
    /// 1D Lagrange shape functions. The
    /// value of all the shape function derivatives at the local coordinate s
    /// are returned in the array DPsi.
    template<unsigned NNODE_1D>
    void d2shape(const double& s, double* DPsi)
    {
      std::ostringstream error_stream;
      error_stream << "One dimensional Lagrange shape function "
                   << "second derivatives "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// 1D shape functions specialised to linear order (2 Nodes)
    // Note that the numbering is such that shape[0] is at s = -1.0.
    // and shape[1] is at s = 1.0
    template<>
    inline void shape<2>(const double& s, double* Psi)
    {
      Psi[0] = 0.5 * (1.0 - s);
      Psi[1] = 0.5 * (1.0 + s);
    }
 
    /// Derivatives of 1D shape functions specialised to linear order (2 Nodes)
    template<>
    inline void dshape<2>(const double& s, double* DPsi)
    {
      DPsi[0] = -0.5;
      DPsi[1] = 0.5;
    }
 
    /// Second Derivatives of 1D shape functions,
    /// specialised to linear order  (2 Nodes)
    template<>
    inline void d2shape<2>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
    }
 
    /// 1D shape functions specialised to quadratic order (3 Nodes)
    // Note that the numbering is such that shape[0] is at s = -1.0,
    // shape[1] is at s = 0.0 and shape[2] is at s = 1.0.
    template<>
    inline void shape<3>(const double& s, double* Psi)
    {
      Psi[0] = 0.5 * s * (s - 1.0);
      Psi[1] = 1.0 - s * s;
      Psi[2] = 0.5 * s * (s + 1.0);
    }
 
    /// Derivatives of 1D shape functions specialised to quadratic order (3
    /// Nodes)
    template<>
    inline void dshape<3>(const double& s, double* DPsi)
    {
      DPsi[0] = s - 0.5;
      DPsi[1] = -2.0 * s;
      DPsi[2] = s + 0.5;
    }
 
 
    /// Second Derivatives of 1D shape functions, specialised to quadratic order
    /// (3 Nodes)
    template<>
    inline void d2shape<3>(const double& s, double* DPsi)
    {
      DPsi[0] = 1.0;
      DPsi[1] = -2.0;
      DPsi[2] = 1.0;
    }
 
    /// 1D shape functions specialised to cubic order (4 Nodes)
    template<>
    inline void shape<4>(const double& s, double* Psi)
    {
      // Output from Maple
      double t1 = s * s;
      double t2 = t1 * s;
      double t3 = 0.5625 * t2;
      double t4 = 0.5625 * t1;
      double t5 = 0.625E-1 * s;
      double t7 = 0.16875E1 * t2;
      double t8 = 0.16875E1 * s;
      Psi[0] = -t3 + t4 + t5 - 0.625E-1;
      Psi[1] = t7 - t4 - t8 + 0.5625;
      Psi[2] = -t7 - t4 + t8 + 0.5625;
      Psi[3] = t3 + t4 - t5 - 0.625E-1;
    }
 
 
    /// Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
    template<>
    inline void dshape<4>(const double& s, double* DPsi)
    {
      // Output from Maple
      double t1 = s * s;
      double t2 = 0.16875E1 * t1;
      double t3 = 0.1125E1 * s;
      double t5 = 0.50625E1 * t1;
      DPsi[0] = -t2 + t3 + 0.625E-1;
      DPsi[1] = t5 - t3 - 0.16875E1;
      DPsi[2] = -t5 - t3 + 0.16875E1;
      DPsi[3] = t2 + t3 - 0.625E-1;
    }
 
    /// Second Derivatives of 1D shape functions specialised to cubic
    /// order (4 Nodes)
    template<>
    inline void d2shape<4>(const double& s, double* DPsi)
    {
      // Output from Maple (modified by ALH, CHECK IT)
      double t1 = 2.0 * s;
      double t2 = 0.16875E1 * t1;
      double t5 = 0.50625E1 * t1;
      DPsi[0] = -t2 + 0.1125E1;
      DPsi[1] = t5 - 0.1125E1;
      DPsi[2] = -t5 - 0.1125E1;
      DPsi[3] = t2 + 0.1125E1;
    }
 
  }; // namespace OneDimLagrange
 
  //======================================================================
  /// One dimensional shape functions and derivatives. Empty -- simply
  /// establishes the template parameters.
  //======================================================================
  namespace OneDimDiscontinuousGalerkin
  {
    /// Definition for 1D Lagrange shape functions. The
    /// value of all the shape functions at the local coordinate s
    /// are returned in the array Psi.
    template<unsigned NNODE_1D>
    void shape(const double& s, double* Psi)
    {
      // Create an output stream
      std::ostringstream error_stream;
 
      // Create the error message
      error_stream << "One dimensional Lagrange shape functions "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
 
      // Throw the error message
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Definition for derivatives of 1D Lagrange shape functions. The
    /// value of all the shape function derivatives at the local coordinate s
    /// are returned in the array DPsi.
    template<unsigned NNODE_1D>
    void dshape(const double& s, double* DPsi)
    {
      // Create an output stream
      std::ostringstream error_stream;
 
      // Create the error message
      error_stream << "One dimensional Lagrange shape function derivatives "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
 
      // Throw the error message
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Definition for second derivatives of
    /// 1D Lagrange shape functions. The
    /// value of all the shape function derivatives at the local coordinate s
    /// are returned in the array DPsi.
    template<unsigned NNODE_1D>
    void d2shape(const double& s, double* DPsi)
    {
      // Create an output stream
      std::ostringstream error_stream;
 
      // Create the error message
      error_stream << "One dimensional Lagrange shape function "
                   << "second derivatives "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
 
      // Throw the error message
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// 1D shape functions specialised to linear order (2 Nodes)
    // Note that the numbering is such that shape[0] is at s = -1.0.
    // and shape[1] is at s = 1.0
    template<>
    inline void shape<2>(const double& s, double* Psi)
    {
      Psi[0] = 0.0;
      Psi[1] = 1.0;
    }
 
    /// Derivatives of 1D shape functions specialised to linear order (2 Nodes)
    template<>
    inline void dshape<2>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
    }
 
    /// Second Derivatives of 1D shape functions,
    /// specialised to linear order  (2 Nodes)
    template<>
    inline void d2shape<2>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
    }
 
    /// 1D shape functions specialised to quadratic order (3 Nodes)
    // Note that the numbering is such that shape[0] is at s = -1.0,
    // shape[1] is at s = 0.0 and shape[2] is at s = 1.0.
    template<>
    inline void shape<3>(const double& s, double* Psi)
    {
      Psi[0] = 0.0;
      Psi[1] = 0.5 * (1.0 - s);
      Psi[2] = 0.5 * (1.0 + s);
    }
 
    /// Derivatives of 1D shape functions specialised to quadratic order (3
    /// Nodes)
    template<>
    inline void dshape<3>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = -0.5;
      DPsi[2] = 0.5;
    }
 
    /// Second Derivatives of 1D shape functions, specialised to quadratic order
    /// (3 Nodes)
    template<>
    inline void d2shape<3>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
      DPsi[2] = 0.0;
    }
 
    /// 1D shape functions specialised to cubic order (4 Nodes)
    template<>
    inline void shape<4>(const double& s, double* Psi)
    {
      Psi[0] = 0.0;
      Psi[1] = 0.5 * s * (s - 1.0);
      Psi[2] = 1.0 - s * s;
      Psi[3] = 0.5 * s * (s + 1.0);
    }
 
 
    /// Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
    template<>
    inline void dshape<4>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = s - 0.5;
      DPsi[2] = -2.0 * s;
      DPsi[3] = s + 0.5;
    }
 
    /// Second derivatives of 1D shape functions specialised to cubic
    /// order (4 nodes)
    template<>
    inline void d2shape<4>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 1.0;
      DPsi[2] = -2.0;
      DPsi[3] = 1.0;
    }
  }; // namespace OneDimDiscontinuousGalerkin
 
 
  //======================================================================
  /// One dimensional shape functions and derivatives. Empty -- simply
  /// establishes the template parameters.
  //======================================================================
  namespace OneDimDiscontinuousGalerkinMixedOrderBasis
  {
    /// Definition for 1D Lagrange shape functions. The
    /// value of all the shape functions at the local coordinate s
    /// are returned in the array Psi.
    template<unsigned NNODE_1D>
    void shape(const double& s, double* Psi)
    {
      // Create an output stream
      std::ostringstream error_stream;
 
      // Create the error message
      error_stream << "One dimensional Lagrange shape functions "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
 
      // Throw the error message
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Definition for derivatives of 1D Lagrange shape functions. The
    /// value of all the shape function derivatives at the local coordinate s
    /// are returned in the array DPsi.
    template<unsigned NNODE_1D>
    void dshape(const double& s, double* DPsi)
    {
      // Create an output stream
      std::ostringstream error_stream;
 
      // Create the error message
      error_stream << "One dimensional Lagrange shape function derivatives "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
 
      // Throw the error message
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Definition for second derivatives of
    /// 1D Lagrange shape functions. The
    /// value of all the shape function derivatives at the local coordinate s
    /// are returned in the array DPsi.
    template<unsigned NNODE_1D>
    void d2shape(const double& s, double* DPsi)
    {
      // Create an output stream
      std::ostringstream error_stream;
 
      // Create the error message
      error_stream << "One dimensional Lagrange shape function "
                   << "second derivatives "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
 
      // Throw the error message
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// 1D shape functions specialised to linear order (2 Nodes)
    // Note that the numbering is such that shape[0] is at s = -1.0.
    // and shape[1] is at s = 1.0
    template<>
    inline void shape<2>(const double& s, double* Psi)
    {
      Psi[0] = 0.5 * (1.0 - s);
      Psi[1] = 0.5 * (1.0 + s);
    }
 
    /// Derivatives of 1D shape functions specialised to linear order (2 Nodes)
    template<>
    inline void dshape<2>(const double& s, double* DPsi)
    {
      DPsi[0] = -0.5;
      DPsi[1] = 0.5;
    }
 
    /// Second Derivatives of 1D shape functions,
    /// specialised to linear order  (2 Nodes)
    template<>
    inline void d2shape<2>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
    }
 
    /// 1D shape functions specialised to quadratic order (3 Nodes)
    // Note that the numbering is such that shape[0] is at s = -1.0,
    // shape[1] is at s = 0.0 and shape[2] is at s = 1.0.
    template<>
    inline void shape<3>(const double& s, double* Psi)
    {
      Psi[0] = 0.5 * (1.0 - s);
      Psi[1] = 0.0;
      Psi[2] = 0.5 * (1.0 + s);
    }
 
    /// Derivatives of 1D shape functions specialised to quadratic order (3
    /// Nodes)
    template<>
    inline void dshape<3>(const double& s, double* DPsi)
    {
      DPsi[0] = -0.5;
      DPsi[1] = 0.0;
      DPsi[2] = 0.5;
    }
 
    /// Second Derivatives of 1D shape functions, specialised to quadratic order
    /// (3 Nodes)
    template<>
    inline void d2shape<3>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
      DPsi[2] = 0.0;
    }
 
    /// 1D shape functions specialised to cubic order (4 Nodes)
    template<>
    inline void shape<4>(const double& s, double* Psi)
    {
      Psi[0] = 0.5 * (1.0 - s);
      Psi[1] = 0.0;
      Psi[2] = 0.0;
      Psi[3] = 0.5 * (1.0 + s);
    }
 
 
    /// Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
    template<>
    inline void dshape<4>(const double& s, double* DPsi)
    {
      DPsi[0] = -0.5;
      DPsi[1] = 0.0;
      DPsi[2] = 0.0;
      DPsi[3] = 0.5;
    }
 
    /// Second derivatives of 1D shape functions specialised to cubic
    /// order (4 nodes)
    template<>
    inline void d2shape<4>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
      DPsi[2] = 0.0;
      DPsi[3] = 0.0;
    }
  }; // namespace OneDimDiscontinuousGalerkinMixedOrderBasis
 
 
  //======================================================================
  /// One dimensional shape functions and derivatives. Empty -- simply
  /// establishes the template parameters.
  //======================================================================
  namespace OneDimDiscontinuousGalerkinMixedOrderTest
  {
    /// Definition for 1D Lagrange shape functions. The
    /// value of all the shape functions at the local coordinate s
    /// are returned in the array Psi.
    template<unsigned NNODE_1D>
    void shape(const double& s, double* Psi)
    {
      // Create an output stream
      std::ostringstream error_stream;
 
      // Create the error message
      error_stream << "One dimensional Lagrange shape functions "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
 
      // Throw the error message
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Definition for derivatives of 1D Lagrange shape functions. The
    /// value of all the shape function derivatives at the local coordinate s
    /// are returned in the array DPsi.
    template<unsigned NNODE_1D>
    void dshape(const double& s, double* DPsi)
    {
      // Create an output stream
      std::ostringstream error_stream;
 
      // Create the error message
      error_stream << "One dimensional Lagrange shape function derivatives "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
 
      // Throw the error message
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// Definition for second derivatives of
    /// 1D Lagrange shape functions. The
    /// value of all the shape function derivatives at the local coordinate s
    /// are returned in the array DPsi.
    template<unsigned NNODE_1D>
    void d2shape(const double& s, double* DPsi)
    {
      // Create an output stream
      std::ostringstream error_stream;
 
      // Create the error message
      error_stream << "One dimensional Lagrange shape function "
                   << "second derivatives "
                   << "have not been defined "
                   << "for " << NNODE_1D << " nodes." << std::endl;
 
      // Throw the error message
      throw OomphLibError(
        error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
    }
 
    /// 1D shape functions specialised to linear order (2 Nodes)
    // Note that the numbering is such that shape[0] is at s = -1.0.
    // and shape[1] is at s = 1.0
    template<>
    inline void shape<2>(const double& s, double* Psi)
    {
      Psi[0] = 0.0;
      Psi[1] = 1.0;
    }
 
    /// Derivatives of 1D shape functions specialised to linear order (2 Nodes)
    template<>
    inline void dshape<2>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
    }
 
    /// Second Derivatives of 1D shape functions,
    /// specialised to linear order  (2 Nodes)
    template<>
    inline void d2shape<2>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
    }
 
    /// 1D shape functions specialised to quadratic order (3 Nodes)
    // Note that the numbering is such that shape[0] is at s = -1.0,
    // shape[1] is at s = 0.0 and shape[2] is at s = 1.0.
    template<>
    inline void shape<3>(const double& s, double* Psi)
    {
      Psi[0] = 0.0;
      Psi[1] = 0.0;
      Psi[2] = 1.0;
    }
 
    /// Derivatives of 1D shape functions specialised to quadratic order (3
    /// Nodes)
    template<>
    inline void dshape<3>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
      DPsi[2] = 0.0;
    }
 
    /// Second Derivatives of 1D shape functions, specialised to quadratic order
    /// (3 Nodes)
    template<>
    inline void d2shape<3>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
      DPsi[2] = 0.0;
    }
 
    /// 1D shape functions specialised to cubic order (4 Nodes)
    template<>
    inline void shape<4>(const double& s, double* Psi)
    {
      Psi[0] = 0.0;
      Psi[1] = 0.0;
      Psi[2] = 0.0;
      Psi[3] = 1.0;
    }
 
 
    /// Derivatives of 1D shape functions specialised to cubic order (4 Nodes)
    template<>
    inline void dshape<4>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
      DPsi[2] = 0.0;
      DPsi[3] = 0.0;
    }
 
    /// Second derivatives of 1D shape functions specialised to cubic
    /// order (4 nodes)
    template<>
    inline void d2shape<4>(const double& s, double* DPsi)
    {
      DPsi[0] = 0.0;
      DPsi[1] = 0.0;
      DPsi[2] = 0.0;
      DPsi[3] = 0.0;
    }
  }; // namespace OneDimDiscontinuousGalerkinMixedOrderTest
 
  //===============================================================
  /// One Dimensional Hermite shape functions
  //===============================================================
  namespace OneDimHermite
  {
    // Convention for polynomial numbering scheme
    // Type 0 is position, 1 is slope
    // Node 0 is at s=0 and 1 is s=1
 
    /// Constructor sets the values of the shape functions at the position s.
    inline void shape(const double& s, double Psi[2][2])
    {
      // Node 0
      Psi[0][0] = 0.25 * (s * s * s - 3.0 * s + 2.0);
      Psi[0][1] = 0.25 * (s * s * s - s * s - s + 1.0);
      // Node 1
      Psi[1][0] = 0.25 * (2.0 + 3.0 * s - s * s * s);
      Psi[1][1] = 0.25 * (s * s * s + s * s - s - 1.0);
    }
 
 
    /// Derivatives of 1D Hermite shape functions
    inline void dshape(const double& s, double DPsi[2][2])
    {
      // Node 0
      DPsi[0][0] = 0.75 * (s * s - 1.0);
      DPsi[0][1] = 0.25 * (3.0 * s * s - 2.0 * s - 1.0);
      // Node 1
      DPsi[1][0] = 0.75 * (1.0 - s * s);
      DPsi[1][1] = 0.25 * (3.0 * s * s + 2.0 * s - 1.0);
    }
 
    /// Second derivatives of the Hermite shape functions
    inline void d2shape(const double& s, double DPsi[2][2])
    {
      // Node 0
      DPsi[0][0] = 1.5 * s;
      DPsi[0][1] = 0.5 * (3.0 * s - 1.0);
      // Node 1
      DPsi[1][0] = -1.5 * s;
      DPsi[1][1] = 0.5 * (3.0 * s + 1.0);
    }
 
  }; // namespace OneDimHermite
 
  //=====================================================================
  /// Class that returns the shape functions associated with legendre
  //=====================================================================
  template<unsigned NNODE_1D>
  class OneDimensionalLegendreShape : public Shape
  {
    static bool Nodes_calculated;
 
  public:
    static Vector<double> z;
 
    /// Static function used to populate the stored positions
    static inline void calculate_nodal_positions()
    {
      if (!Nodes_calculated)
      {
        Orthpoly::gll_nodes(NNODE_1D, z);
        Nodes_calculated = true;
      }
    }
 
    static inline double nodal_position(const unsigned& n)
    {
      return z[n];
    }
 
    /// Constructor
    OneDimensionalLegendreShape(const double& s) : Shape(NNODE_1D)
    {
      using namespace Orthpoly;
 
      unsigned p = NNODE_1D - 1;
      // Now populate the shape function
      for (unsigned i = 0; i < NNODE_1D; i++)
      {
        // If we're at one of the nodes, the value must be 1.0
        if (std::fabs(s - z[i]) < Orthpoly::eps)
        {
          (*this)[i] = 1.0;
        }
        // Otherwise use the lagrangian interpolant
        else
        {
          (*this)[i] = (1.0 - s * s) * dlegendre(p, s) /
                       (p * (p + 1) * legendre(p, z[i]) * (z[i] - s));
        }
      }
    }
  };
 
  template<unsigned NNODE_1D>
  Vector<double> OneDimensionalLegendreShape<NNODE_1D>::z;
 
  template<unsigned NNODE_1D>
  bool OneDimensionalLegendreShape<NNODE_1D>::Nodes_calculated = false;
 
 
  template<unsigned NNODE_1D>
  class OneDimensionalLegendreDShape : public Shape
  {
  public:
    // Constructor
    OneDimensionalLegendreDShape(const double& s) : Shape(NNODE_1D)
    {
      unsigned p = NNODE_1D - 1;
      Vector<double> z = OneDimensionalLegendreShape<NNODE_1D>::z;
 
 
      bool root = false;
 
      for (unsigned i = 0; i < NNODE_1D; i++)
      {
        unsigned rootnum = 0;
        for (unsigned j = 0; j < NNODE_1D; j++)
        { // Loop over roots to check if
          if (std::fabs(s - z[j]) < 10 * Orthpoly::eps)
          { // s happens to be a root.
            root = true;
            break;
          }
          rootnum += 1;
        }
        if (root == true)
        {
          if (i == rootnum && i == 0)
          {
            (*this)[i] = -(1.0 + p) * p / 4.0;
          }
          else if (i == rootnum && i == p)
          {
            (*this)[i] = (1.0 + p) * p / 4.0;
          }
          else if (i == rootnum)
          {
            (*this)[i] = 0.0;
          }
          else
          {
            (*this)[i] = Orthpoly::legendre(p, z[rootnum]) /
                         Orthpoly::legendre(p, z[i]) / (z[rootnum] - z[i]);
          }
        }
        else
        {
          (*this)[i] =
            ((1 + s * (s - 2 * z[i])) / (s - z[i]) * Orthpoly::dlegendre(p, s) -
             (1 - s * s) * Orthpoly::ddlegendre(p, s)) /
            p / (p + 1.0) / Orthpoly::legendre(p, z[i]) / (s - z[i]);
        }
        root = false;
      }
    }
  };
 
 
  //=====================================================================
  /// Non-templated class that returns modal hierachical shape functions
  /// based on Legendre polynomials
  //=====================================================================
  class OneDimensionalModalShape : public Shape
  {
  public:
    /// Constructor
    OneDimensionalModalShape(const unsigned p_order, const double& s)
      : Shape(p_order)
    {
      // Populate the shape functions
      (*this)[0] = 0.5 * (1.0 - s);
      (*this)[1] = 0.5 * (1.0 + s);
      for (unsigned i = 2; i < p_order; i++)
      {
        (*this)[i] =
          (0.5 * (1.0 - s)) * (0.5 * (1.0 + s)) * Orthpoly::legendre(i - 2, s);
      }
    }
  };
 
  class OneDimensionalModalDShape : public Shape
  {
  public:
    // Constructor
    OneDimensionalModalDShape(const unsigned p_order, const double& s)
      : Shape(p_order)
    {
      // Populate the shape functions
      (*this)[0] = -0.5;
      (*this)[1] = 0.5;
      for (unsigned i = 2; i < p_order; i++)
      {
        (*this)[i] = (0.5 * (1.0 - s)) * (0.5 * (1.0 + s)) *
                       Orthpoly::dlegendre(i - 2, s) -
                     0.5 * s * Orthpoly::legendre(i - 2, s);
      }
    }
  };
 
} // namespace oomph
 
#endif