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-2026 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
#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;
 
    /// Boolean to indicate whether the data has been copied
    bool Is_psi_copied;
 
    /// Boolean to indicate whether the data is a copy
    bool Is_psi_a_copy;
 
    /// 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), Is_psi_copied(false), Is_psi_a_copy(false)
    {
      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), Is_psi_copied(false), Is_psi_a_copy(false)
    {
      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),
        Is_psi_copied(false),
        Is_psi_a_copy(false)
    {
    }
 
    /// Assignment operator, shallow copy (not recommended so warning given)
    void operator=(Shape& shape)
    {
      if (!Suppress_warning_about_assignment)
      {
        std::ostringstream warning_stream;
        warning_stream << "The assignment operator, =, makes a shallow copy."
                       << std::endl
                       << "This can have unexpected consequences, so we"
                       << std::endl
                       << "now recommend using the explicitly named function "
                       << std::endl
                       << "Shape::shallow_copy_from(.)." << std::endl
                       << "If you want to suppress" << std::endl
                       << "this warning then set the static boolean:"
                       << std::endl
                       << "Shape::Suppress_warning_about_assignment"
                       << std::endl
                       << " to true." << std::endl;
 
        OomphLibWarning(warning_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
 
      this->shallow_copy_from(shape);
    }
 
    /// Assignment operator, shallow copy (not recommended so warning given)
    void operator=(Shape* const& shape_pt)
    {
      if (!Suppress_warning_about_assignment)
      {
        std::ostringstream warning_stream;
        warning_stream << "The assignment operator, =, makes a shallow copy."
                       << std::endl
                       << "This can have unexpected consequences, so we"
                       << std::endl
                       << "now recommend using the explicitly named function "
                       << std::endl
                       << "Shape::shallow_copy_from(.)." << std::endl
                       << "If you want to suppress" << std::endl
                       << "this warning then set the static boolean:"
                       << std::endl
                       << "Shape::Suppress_warning_about_assignment"
                       << std::endl
                       << " to true." << std::endl;
 
        OomphLibWarning(warning_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
 
      this->shallow_copy_from(shape_pt);
    }
 
    /// This function makes a shallow copy
    /// (sets the pointer to the allocated data to be from another Shape object)
    void shallow_copy_from(Shape& shape)
    {
      // Sets the indices of this object to be consistent with that
      // being copied
      Index1 = shape.nindex1();
      Index2 = shape.nindex2();
      // Set the pointer to the data to be that from the copied object
      Psi = shape.Psi;
      // Indicate that data has been copied
      shape.Is_psi_copied = true;
      Is_psi_a_copy = true;
      // Delete own storage??
    }
 
 
    /// This function make a shallow copy
    /// (sets the pointer to the allocated data to be from another Shape object)
    void shallow_copy_from(Shape* const& shape_pt)
    {
      // Sets the indices of this object to be consistent with that
      // being copied
      Index1 = shape_pt->nindex1();
      Index2 = shape_pt->nindex2();
      // Set pointer
      Psi = shape_pt->Psi;
      // Indicate that data has been copied
      shape_pt->Is_psi_copied = true;
      Is_psi_a_copy = true;
    }
 
    /// Destructor, clear up the memory allocated by the object
    ~Shape()
    {
      delete[] Allocated_storage;
      Allocated_storage = 0;
    }
 
  private:
    /// Change the size of the storage
    void resize(const unsigned& N, const unsigned& M = 1)
    {
#ifdef PARANOID
      // If the storage is copied, do not allow a resize
      if (Is_psi_copied)
      {
        std::ostringstream error_stream;
        error_stream << "Cannot resize an object whose storage has been copied."
                     << std::endl
                     << "This could lead to a dangling pointer in the copy."
                     << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
      // 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;
    }
 
  public:
    /// 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;
    }
 
    /// Boolean used to suppress warning about assignment operator
    static bool Suppress_warning_about_assignment;
  };
 
  //================================================================
  /// 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;
 
    /// Boolean to indicate whether the data has been copied
    bool Is_dpsi_copied;
 
    /// Boolean to indicate whether the data is a copy
    bool Is_dpsi_a_copy;
 
 
    /// 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),
        Is_dpsi_copied(false),
        Is_dpsi_a_copy(false)
    {
      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),
        Is_dpsi_copied(false),
        Is_dpsi_a_copy(false)
    {
      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),
        Is_dpsi_copied(false),
        Is_dpsi_a_copy(false)
    {
    }
 
    /// Broken copy constructor
    DShape(const DShape& dshape) = delete;
 
    /// Assignment operator, shallow copy (not recommended so warning given)
    void operator=(DShape& dshape)
    {
      if (!Suppress_warning_about_assignment)
      {
        std::ostringstream warning_stream;
        warning_stream << "The assignment operator, =, makes a shallow copy."
                       << std::endl
                       << "This can have unexpected consequences, so we"
                       << std::endl
                       << "now recommend using the explicitly named function "
                       << std::endl
                       << "DShape::shallow_copy_from(.)." << std::endl
                       << "If you want to suppress" << std::endl
                       << "this warning then set the static boolean:"
                       << std::endl
                       << "DShape::Suppress_warning_about_assignment"
                       << std::endl
                       << " to true." << std::endl;
 
        OomphLibWarning(warning_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
 
      this->shallow_copy_from(dshape);
    }
 
 
    /// Assignment operator, shallow copy (not reccomended so warning given)
    void operator=(DShape* const& dshape_pt)
    {
      if (!Suppress_warning_about_assignment)
      {
        std::ostringstream warning_stream;
        warning_stream << "The assignment operator, =, makes a shallow copy."
                       << std::endl
                       << "This can have unexpected consequences, so we"
                       << std::endl
                       << "now recommend using the explicitly named function "
                       << std::endl
                       << "DShape::shallow_copy_from(.)." << std::endl
                       << "If you want to suppress" << std::endl
                       << "this warning then set the static boolean:"
                       << std::endl
                       << "DShape::Suppress_warning_about_assignment"
                       << std::endl
                       << " to true." << std::endl;
 
        OomphLibWarning(warning_stream.str(),
                        OOMPH_CURRENT_FUNCTION,
                        OOMPH_EXCEPTION_LOCATION);
      }
 
      this->shallow_copy_from(dshape_pt);
    }
 
 
    /// This function does a shallow copy
    /// (resets the pointer to the data)
    void shallow_copy_from(DShape& dshape)
    {
      // Set the values of the indices from the source object
      Index1 = dshape.nindex1();
      Index2 = dshape.nindex2();
      Index3 = dshape.nindex3();
      // Copy the pointer and mark is as being a copy
      DPsi = dshape.DPsi;
      Is_dpsi_a_copy = true;
      dshape.Is_dpsi_copied = true;
    }
 
    /// This function does a shallow copy
    /// (resets the pointer to the data)
    void shallow_copy_from(DShape* const& dshape_pt)
    {
      // Set values of the indices from the source object
      Index1 = dshape_pt->nindex1();
      Index2 = dshape_pt->nindex2();
      Index3 = dshape_pt->nindex3();
      // Copy pointer and mark object as being a copy
      DPsi = dshape_pt->DPsi;
      Is_dpsi_a_copy = true;
      dshape_pt->Is_dpsi_copied = true;
    }
 
 
    /// Destructor, clean up the memory allocated by this object
    ~DShape()
    {
      delete[] Allocated_storage;
      Allocated_storage = 0;
    }
 
  private:
    /// 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)
    {
#ifdef PARANOID
      // If the storage is copied, do not allow a resize
      if (Is_dpsi_copied)
      {
        std::ostringstream error_stream;
        error_stream << "Cannot resize an object whose storage has been copied."
                     << std::endl
                     << "This could lead to a dangling pointer in the copy."
                     << std::endl;
        throw OomphLibError(
          error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
      }
#endif
 
      // 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;
    }
 
  public:
    /// 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;
    }
 
    /// Boolean used to suppress warning about assignment operator
    static bool Suppress_warning_about_assignment;
  };
 
  ////////////////////////////////////////////////////////////////////
  //
  // 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