integral.h

Go to the documentation of this file.

// LIC// ====================================================================
// LIC// This file forms part of oomph-lib, the object-oriented,
// LIC// multi-physics finite-element library, available
// LIC// at http://www.oomph-lib.org.
// LIC//
// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
// LIC//
// LIC// This library is free software; you can redistribute it and/or
// LIC// modify it under the terms of the GNU Lesser General Public
// LIC// License as published by the Free Software Foundation; either
// LIC// version 2.1 of the License, or (at your option) any later version.
// LIC//
// LIC// This library is distributed in the hope that it will be useful,
// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// LIC// Lesser General Public License for more details.
// LIC//
// LIC// You should have received a copy of the GNU Lesser General Public
// LIC// License along with this library; if not, write to the Free Software
// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
// LIC// 02110-1301  USA.
// LIC//
// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
// LIC//
// LIC//====================================================================
// Header file for numerical integration routines based on quadrature
 
// Include guards to prevent multiple inclusions of the header
#ifndef OOMPH_INTEGRAL_HEADER
#define OOMPH_INTEGRAL_HEADER
 
 
// Config header generated by autoconfig
#ifdef HAVE_CONFIG_H
#include <oomph-lib-config.h>
#endif
 
// oomph-lib headers
#include "oomph_utilities.h"
#include "orthpoly.h"
 
namespace oomph
{
  //=========================================================
  //=========================================================
  class Integral
  {
  public:
    Integral(){};
 
    Integral(const Integral& dummy) = delete;
 
    void operator=(const Integral&) = delete;
 
    virtual ~Integral(){};
 
    virtual unsigned nweight() const = 0;
 
    virtual double knot(const unsigned& i, const unsigned& j) const = 0;
 
    virtual Vector<double> knot(const unsigned& i) const
    {
      throw OomphLibError("Not implemented for this integration scheme (yet?).",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
    }
 
    virtual double weight(const unsigned& i) const = 0;
  };
 
  //=============================================================================
  //=============================================================================
  class PointIntegral : public Integral
  {
  public:
    PointIntegral(){};
 
    PointIntegral(const PointIntegral& dummy) = delete;
 
    void operator=(const PointIntegral&) = delete;
 
    unsigned nweight() const
    {
      return 1;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      throw OomphLibError("Local coordinate vector is of size zero, so this "
                          "should never be called.",
                          OOMPH_CURRENT_FUNCTION,
                          OOMPH_EXCEPTION_LOCATION);
 
      // Dummy return
      return 0.0;
    }
 
    double weight(const unsigned& i) const
    {
      return 1.0;
    }
  };
 
  //=========================================================
  //=========================================================
  template<unsigned DIM, unsigned NPTS_1D>
  class Gauss
  {
  };
 
 
  //=========================================================
  //=========================================================
  template<>
  class Gauss<1, 2> : public Integral
  {
  private:
    static const unsigned Npts = 2;
    static const double Knot[2][1], Weight[2];
 
  public:
    Gauss(){};
 
    Gauss(const Gauss& dummy) = delete;
 
    void operator=(const Gauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class Gauss<1, 3> : public Integral
  {
  private:
    static const unsigned Npts = 3;
    static const double Knot[3][1], Weight[3];
 
  public:
    Gauss(){};
 
    Gauss(const Gauss& dummy) = delete;
 
    void operator=(const Gauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class Gauss<1, 4> : public Integral
  {
  private:
    static const unsigned Npts = 4;
    static const double Knot[4][1], Weight[4];
 
  public:
    Gauss(){};
 
    Gauss(const Gauss& dummy) = delete;
 
    void operator=(const Gauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class Gauss<2, 2> : public Integral
  {
  private:
    static const unsigned Npts = 4;
    static const double Knot[4][2], Weight[4];
 
  public:
    Gauss(){};
 
    Gauss(const Gauss& dummy) = delete;
 
    void operator=(const Gauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class Gauss<2, 3> : public Integral
  {
  private:
    static const unsigned Npts = 9;
    static const double Knot[9][2], Weight[9];
 
  public:
    Gauss(){};
 
    Gauss(const Gauss& dummy) = delete;
 
    void operator=(const Gauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class Gauss<2, 4> : public Integral
  {
  private:
    static const unsigned Npts = 16;
    static const double Knot[16][2], Weight[16];
 
 
  public:
    Gauss(){};
 
    Gauss(const Gauss& dummy) = delete;
 
    void operator=(const Gauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class Gauss<3, 2> : public Integral
  {
  private:
    static const unsigned Npts = 8;
    static const double Knot[8][3], Weight[8];
 
  public:
    Gauss(){};
 
    Gauss(const Gauss& dummy) = delete;
 
    void operator=(const Gauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class Gauss<3, 3> : public Integral
  {
  private:
    static const unsigned Npts = 27;
    static const double Knot[27][3], Weight[27];
 
  public:
    Gauss(){};
 
    Gauss(const Gauss& dummy) = delete;
 
    void operator=(const Gauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class Gauss<3, 4> : public Integral
  {
  private:
    static const unsigned Npts = 64;
    static const double Knot[64][3], Weight[64];
 
  public:
    Gauss(){};
 
    Gauss(const Gauss& dummy) = delete;
 
    void operator=(const Gauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<unsigned DIM, unsigned NPTS_1D>
  class Gauss_Rescaled : public Gauss<DIM, NPTS_1D>
  {
  private:
    double Lower, Upper, Range;
 
  public:
    Gauss_Rescaled(){};
 
    Gauss_Rescaled(const Gauss_Rescaled& dummy) = delete;
 
    void operator=(const Gauss_Rescaled&) = delete;
 
    Gauss_Rescaled(double lower, double upper) : Lower(lower), Upper(upper)
    {
      // Set the range of integration
      Range = upper - lower;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return (0.5 * (Gauss<DIM, NPTS_1D>::knot(i, j) * Range + Lower + Upper));
    }
 
    double weight(const unsigned& i) const
    {
      return Gauss<DIM, NPTS_1D>::weight(i) *
             pow(0.5 * Range, static_cast<int>(DIM));
    }
  };
 
  //=========================================================
  //=========================================================
  template<unsigned DIM, unsigned NPTS_1D>
  class TGauss
  {
  };
 
 
  //=========================================================
  //=========================================================
  template<>
  class TGauss<1, 2> : public Integral
  {
  private:
    static const unsigned Npts = 2;
 
    static const double Knot[2][1], Weight[2];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class TGauss<1, 3> : public Integral
  {
  private:
    static const unsigned Npts = 3;
    static const double Knot[3][1], Weight[3];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
 
  //=========================================================
  //=========================================================
  template<>
  class TGauss<1, 4> : public Integral
  {
  private:
    static const unsigned Npts = 4;
    static const double Knot[4][1], Weight[4];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  template<>
  class TGauss<1, 5> : public Integral
  {
  private:
    static const unsigned Npts = 5;
    static const double Knot[5][1], Weight[5];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
 
  //=========================================================
  //=========================================================
  template<>
  class TGauss<2, 2> : public Integral
  {
  private:
    static const unsigned Npts = 3;
 
    static const double Knot[3][2], Weight[3];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class TGauss<2, 3> : public Integral
  {
  private:
    static const unsigned Npts = 7;
    static const double Knot[7][2], Weight[7];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
 
  //=========================================================
  //=========================================================
  template<>
  class TGauss<2, 4> : public Integral
  {
  private:
    static const unsigned Npts = 13;
    static const double Knot[13][2], Weight[13];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //------------------------------------------------------------
  //"Full integration" weights for 2D triangles
  // accurate up to order 11
  // http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/
  // quadrature_rules_tri.html
  //------------------------------------------------------------
  template<>
  class TGauss<2, 13> : public Integral
  {
  private:
    static const unsigned Npts = 37;
    static const double Knot[37][2], Weight[37];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //------------------------------------------------------------
  //"Full integration" weights for 2D triangles
  // accurate up to points 19, degree of precision 8, a rule from ACM TOMS
  // algorithm #584.
  // http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/quadrature_rules_tri.html
  // TOMS589_19 a 19 point Gauss rule accurate to order 8
  //------------------------------------------------------------
  template<>
  class TGauss<2, 9> : public Integral
  {
  private:
    static const unsigned Npts = 19;
    static const double Knot[19][2], Weight[19];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //------------------------------------------------------------
  //"Full integration" weights for 2D triangles
  // accurate up to order 16
  // http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/
  // quadrature_rules_tri.html
  // This uses the order 16 Dunavant rule, computed using software available
  // from
  // https://people.sc.fsu.edu/~jburkardt/cpp_src/triangle_dunavant_rule/triangle_dunavant_rule.html
  // This integegration scheme is accurate up to order 16 and uses 52 points
  //------------------------------------------------------------
  template<>
  class TGauss<2, 16> : public Integral
  {
  private:
    static const unsigned Npts = 52;
    static const double Knot[52][2], Weight[52];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //------------------------------------------------------------
  //"Full integration" weights for 2D triangles
  // accurate up to order 15
  // http://people.sc.fsu.edu/~jburkardt/datasets/quadrature_rules_tri/quadrature_rules_tri.html
  //------------------------------------------------------------
 
  template<>
  class TGauss<2, 5> : public Integral
  {
  private:
    static const unsigned Npts = 64;
    static const double Knot[64][2], Weight[64];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=========================================================
  //=========================================================
  template<>
  class TGauss<3, 2> : public Integral
  {
  private:
    static const unsigned Npts = 4;
    static const double Knot[4][3], Weight[4];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
 
  //=========================================================
  //=========================================================
  template<>
  class TGauss<3, 3> : public Integral
  {
  private:
    static const unsigned Npts = 11;
    static const double Knot[11][3], Weight[11];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
 
  //=========================================================
  //=========================================================
  template<>
  class TGauss<3, 5> : public Integral
  {
  private:
    static const unsigned Npts = 45;
    static const double Knot[45][3], Weight[45];
 
  public:
    TGauss(){};
 
    TGauss(const TGauss& dummy) = delete;
 
    void operator=(const TGauss&) = delete;
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
 
  //===================================================================
  //===================================================================
  template<unsigned DIM, unsigned NPTS_1D>
  class GaussLobattoLegendre
  {
  };
 
  //===================================================================
  //===================================================================
  template<unsigned NPTS_1D>
  class GaussLobattoLegendre<1, NPTS_1D> : public Integral
  {
  private:
    static const unsigned Npts = NPTS_1D;
    // These are not constant, because they are calculated on the fly
    double Knot[NPTS_1D][1], Weight[NPTS_1D];
 
  public:
    GaussLobattoLegendre();
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
 
  //=============================================================
  //=============================================================
  template<unsigned NPTS_1D>
  GaussLobattoLegendre<1, NPTS_1D>::GaussLobattoLegendre()
  {
    // Temporary storage for the integration points
    Vector<double> s(NPTS_1D), w(NPTS_1D);
    // call the function that calculates the points
    Orthpoly::gll_nodes(NPTS_1D, s, w);
    // Populate the arrays
    for (unsigned i = 0; i < NPTS_1D; i++)
    {
      Knot[i][0] = s[i];
      Weight[i] = w[i];
    }
  }
 
 
  //===================================================================
  //===================================================================
  template<unsigned NPTS_1D>
  class GaussLobattoLegendre<2, NPTS_1D> : public Integral
  {
  private:
    static const unsigned long int Npts = NPTS_1D * NPTS_1D;
 
    double Knot[NPTS_1D * NPTS_1D][2],
      Weight[NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
    // const BECAUSE THEY ARE CALCULATED (at least currently)
 
  public:
    GaussLobattoLegendre();
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=============================================================
  //=============================================================
 
  template<unsigned NPTS_1D>
  GaussLobattoLegendre<2, NPTS_1D>::GaussLobattoLegendre()
  {
    // Tempoarary storage for the 1D knots and weights
    Vector<double> s(NPTS_1D), w(NPTS_1D);
    // Call the function to populate the array
    Orthpoly::gll_nodes(NPTS_1D, s, w);
    int i_fast = 0, i_slow = 0;
    for (unsigned i = 0; i < NPTS_1D * NPTS_1D; i++)
    {
      if (i_fast == NPTS_1D)
      {
        i_fast = 0;
        i_slow++;
      }
      Knot[i][0] = s[i_fast];
      Knot[i][1] = s[i_slow];
      Weight[i] = w[i_fast] * w[i_slow];
      i_fast++;
    }
  }
 
 
  //===================================================================
  //===================================================================
  template<unsigned NPTS_1D>
  class GaussLobattoLegendre<3, NPTS_1D> : public Integral
  {
  private:
    static const unsigned long int Npts = NPTS_1D * NPTS_1D * NPTS_1D;
 
    double Knot[NPTS_1D * NPTS_1D * NPTS_1D][3],
      Weight[NPTS_1D * NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
    // const BECAUSE THEY ARE CALCULATED (at least currently)
 
  public:
    GaussLobattoLegendre();
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=============================================================
  //=============================================================
 
  template<unsigned NPTS_1D>
  GaussLobattoLegendre<3, NPTS_1D>::GaussLobattoLegendre()
  {
    // Tempoarary storage for the 1D knots and weights
    Vector<double> s(NPTS_1D), w(NPTS_1D);
    // Call the function to populate the array
    Orthpoly::gll_nodes(NPTS_1D, s, w);
    for (unsigned k = 0; k < NPTS_1D; k++)
    {
      for (unsigned j = 0; j < NPTS_1D; j++)
      {
        for (unsigned i = 0; i < NPTS_1D; i++)
        {
          unsigned index = NPTS_1D * NPTS_1D * k + NPTS_1D * j + i;
          Knot[index][0] = s[i];
          Knot[index][1] = s[j];
          Knot[index][2] = s[k];
          Weight[index] = w[i] * w[j] * w[k];
        }
      }
    }
  }
 
 
 
  //===================================================================
  //===================================================================
  template<unsigned DIM, unsigned NPTS_1D>
  class GaussLegendre
  {
  };
 
  //===================================================================
  //===================================================================
  template<unsigned NPTS_1D>
  class GaussLegendre<1, NPTS_1D> : public Integral
  {
  private:
    static const unsigned Npts = NPTS_1D;
    // These are not constant, because they are calculated on the fly
    double Knot[NPTS_1D][1], Weight[NPTS_1D];
 
  public:
    GaussLegendre();
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
 
  //=============================================================
  //=============================================================
  template<unsigned NPTS_1D>
  GaussLegendre<1, NPTS_1D>::GaussLegendre()
  {
    // Temporary storage for the integration points
    Vector<double> s(NPTS_1D), w(NPTS_1D);
    // call the function that calculates the points
    Orthpoly::gl_nodes(NPTS_1D, s, w);
    // Populate the arrays
    for (unsigned i = 0; i < NPTS_1D; i++)
    {
      Knot[i][0] = s[i];
      Weight[i] = w[i];
    }
  }
 
 
  //===================================================================
  //===================================================================
  template<unsigned NPTS_1D>
  class GaussLegendre<2, NPTS_1D> : public Integral
  {
  private:
    static const unsigned long int Npts = NPTS_1D * NPTS_1D;
 
    double Knot[NPTS_1D * NPTS_1D][2],
      Weight[NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
    // const BECAUSE THEY ARE CALCULATED (at least currently)
 
  public:
    GaussLegendre();
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=============================================================
  //=============================================================
 
  template<unsigned NPTS_1D>
  GaussLegendre<2, NPTS_1D>::GaussLegendre()
  {
    // Tempoarary storage for the 1D knots and weights
    Vector<double> s(NPTS_1D), w(NPTS_1D);
    // Call the function to populate the array
    Orthpoly::gl_nodes(NPTS_1D, s, w);
    int i_fast = 0, i_slow = 0;
    for (unsigned i = 0; i < NPTS_1D * NPTS_1D; i++)
    {
      if (i_fast == NPTS_1D)
      {
        i_fast = 0;
        i_slow++;
      }
      Knot[i][0] = s[i_fast];
      Knot[i][1] = s[i_slow];
      Weight[i] = w[i_fast] * w[i_slow];
      i_fast++;
    }
  }
 
 
  //===================================================================
  //===================================================================
  template<unsigned NPTS_1D>
  class GaussLegendre<3, NPTS_1D> : public Integral
  {
  private:
    static const unsigned long int Npts = NPTS_1D * NPTS_1D * NPTS_1D;
 
    double Knot[NPTS_1D * NPTS_1D * NPTS_1D][3],
      Weight[NPTS_1D * NPTS_1D * NPTS_1D]; // COULDN'T MAKE THESE
    // const BECAUSE THEY ARE CALCULATED (at least currently)
 
  public:
    GaussLegendre();
 
    unsigned nweight() const
    {
      return Npts;
    }
 
    double knot(const unsigned& i, const unsigned& j) const
    {
      return Knot[i][j];
    }
 
    double weight(const unsigned& i) const
    {
      return Weight[i];
    }
  };
 
  //=============================================================
  //=============================================================
 
  template<unsigned NPTS_1D>
  GaussLegendre<3, NPTS_1D>::GaussLegendre()
  {
    // Tempoarary storage for the 1D knots and weights
    Vector<double> s(NPTS_1D), w(NPTS_1D);
    // Call the function to populate the array
    Orthpoly::gl_nodes(NPTS_1D, s, w);
    for (unsigned k = 0; k < NPTS_1D; k++)
    {
      for (unsigned j = 0; j < NPTS_1D; j++)
      {
        for (unsigned i = 0; i < NPTS_1D; i++)
        {
          unsigned index = NPTS_1D * NPTS_1D * k + NPTS_1D * j + i;
          Knot[index][0] = s[i];
          Knot[index][1] = s[j];
          Knot[index][2] = s[k];
          Weight[index] = w[i] * w[j] * w[k];
        }
      }
    }
  }
 
 
} // namespace oomph
 
#endif