PolynomialQuadrature.h

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#pragma once

#include "libmesh/utility.h"

#include "DataIO.h"
#include "Uniform.h"
#include "Normal.h"
#include "CartesianProduct.h"
#include "nlohmann/json.h"

namespace PolynomialQuadrature
{
class Polynomial;

std::unique_ptr<const Polynomial> makePolynomial(const Distribution * dist);

class Polynomial
{
public:
  Polynomial() {}
  virtual ~Polynomial() = default;
  virtual void store(std::ostream & stream, void * context) const;
  virtual void store(nlohmann::json & json) const;
  virtual Real compute(const unsigned int order, const Real x, const bool normalize = true) const;
  virtual Real
  computeDerivative(const unsigned int order, const Real x, const unsigned int m = 1) const;
  virtual Real innerProduct(const unsigned int order) const = 0;

  virtual void gaussQuadrature(const unsigned int order,
                               std::vector<Real> & points,
                               std::vector<Real> & weights) const = 0;
  virtual void clenshawQuadrature(const unsigned int order,
                                  std::vector<Real> & points,
                                  std::vector<Real> & weights) const;
  Real productIntegral(const std::vector<unsigned int> order) const;
};

class Legendre : public Polynomial
{
public:
  Legendre(const Real lower_bound, const Real upper_bound);
  virtual void store(std::ostream & stream, void * context) const override;
  virtual void store(nlohmann::json & json) const override;

  virtual Real
  compute(const unsigned int order, const Real x, const bool normalize = true) const override;
  virtual Real computeDerivative(const unsigned int order,
                                 const Real x,
                                 const unsigned int m = 1) const override;
  Real computeDerivativeRef(const unsigned int order, const Real x, const unsigned int m = 1) const;
  virtual Real innerProduct(const unsigned int order) const override;

  virtual void gaussQuadrature(const unsigned int order,
                               std::vector<Real> & points,
                               std::vector<Real> & weights) const override;

  virtual void clenshawQuadrature(const unsigned int order,
                                  std::vector<Real> & points,
                                  std::vector<Real> & weights) const override;

private:
  const Real _lower_bound;
  const Real _upper_bound;
};

class Hermite : public Polynomial
{
public:
  Hermite(const Real mu, const Real sig);
  virtual void store(std::ostream & stream, void * context) const override;
  virtual void store(nlohmann::json & json) const override;

  virtual Real
  compute(const unsigned int order, const Real x, const bool normalize = true) const override;

  virtual Real computeDerivative(const unsigned int order,
                                 const Real x,
                                 const unsigned int m = 1) const override;
  virtual Real innerProduct(const unsigned int order) const override;

  virtual void gaussQuadrature(const unsigned int order,
                               std::vector<Real> & points,
                               std::vector<Real> & weights) const override;

private:
  const Real _mu;
  const Real _sig;
};

class Quadrature
{
public:
  Quadrature() {}
  virtual ~Quadrature() = default;

  virtual unsigned int nPoints() const = 0;
  virtual unsigned int nDim() const = 0;

  virtual std::vector<Real> quadraturePoint(const unsigned int n) const = 0;
  virtual Real quadraturePoint(const unsigned int n, const unsigned int dim) const = 0;
  virtual Real quadratureWeight(const unsigned int n) const = 0;
};

class TensorGrid : public Quadrature
{
public:
  TensorGrid(const std::vector<unsigned int> & npoints,
             std::vector<std::unique_ptr<const Polynomial>> & poly);
  TensorGrid(const unsigned int npoints, std::vector<std::unique_ptr<const Polynomial>> & poly)
    : TensorGrid(std::vector<unsigned int>(poly.size(), npoints), poly)
  {
  }

  virtual unsigned int nPoints() const override { return _quad->numRows(); };
  virtual unsigned int nDim() const override { return _quad->numCols(); };

  virtual std::vector<Real> quadraturePoint(const unsigned int n) const override
  {
    return _quad->computeRow(n);
  };
  virtual Real quadraturePoint(const unsigned int n, const unsigned int dim) const override
  {
    return _quad->computeValue(n, dim);
  };
  virtual Real quadratureWeight(const unsigned int n) const override
  {
    return _quad->computeWeight(n);
  };

private:
  std::unique_ptr<const StochasticTools::WeightedCartesianProduct<Real, Real>> _quad = nullptr;
};

class SmolyakGrid : public Quadrature
{
public:
  SmolyakGrid(const unsigned int max_order, std::vector<std::unique_ptr<const Polynomial>> & poly);

  virtual unsigned int nPoints() const override { return _npoints.back(); }
  virtual unsigned int nDim() const override { return _ndim; }

  virtual std::vector<Real> quadraturePoint(const unsigned int n) const override;
  virtual Real quadraturePoint(const unsigned int n, const unsigned int dim) const override;
  virtual Real quadratureWeight(const unsigned int n) const override;

private:
  unsigned int gridIndex(const unsigned int n) const;
  const unsigned int _ndim;
  std::vector<std::size_t> _npoints;
  std::vector<int> _coeff;
  std::vector<std::unique_ptr<const StochasticTools::WeightedCartesianProduct<Real, Real>>> _quad;
};

class ClenshawCurtisGrid : public Quadrature
{
public:
  ClenshawCurtisGrid(const unsigned int max_order,
                     std::vector<std::unique_ptr<const Polynomial>> & poly);

  virtual unsigned int nPoints() const override { return _quadrature.size(); }
  virtual unsigned int nDim() const override { return _ndim; }

  virtual std::vector<Real> quadraturePoint(const unsigned int n) const override
  {
    return _quadrature[n].first;
  }
  virtual Real quadraturePoint(const unsigned int n, const unsigned int dim) const override
  {
    return _quadrature[n].first[dim];
  }
  virtual Real quadratureWeight(const unsigned int n) const override
  {
    return _quadrature[n].second;
  }

private:
  const unsigned int _ndim;
  std::vector<std::pair<std::vector<Real>, Real>> _quadrature;
};

Real legendre(const unsigned int order,
              const Real x,
              const Real lower_bound = -1.0,
              const Real upper_bound = 1.0);

Real hermite(const unsigned int order, const Real x, const Real mu = 0.0, const Real sig = 1.0);

void gauss_legendre(const unsigned int order,
                    std::vector<Real> & points,
                    std::vector<Real> & weights,
                    const Real lower_bound,
                    const Real upper_bound);

void gauss_hermite(const unsigned int order,
                   std::vector<Real> & points,
                   std::vector<Real> & weights,
                   const Real mu,
                   const Real sig);

void
clenshaw_curtis(const unsigned int order, std::vector<Real> & points, std::vector<Real> & weights);
}

template <>
void dataStore(std::ostream & stream,
               std::unique_ptr<const PolynomialQuadrature::Polynomial> & ptr,
               void * context);
template <>
void dataLoad(std::istream & stream,
              std::unique_ptr<const PolynomialQuadrature::Polynomial> & ptr,
              void * context);