quadrature.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_quadrature_h
#define dealii_quadrature_h


#include <deal.II/base/config.h>
#include <deal.II/base/point.h>
#include <deal.II/base/subscriptor.h>
#include <vector>
#include <array>
#include <memory>

DEAL_II_NAMESPACE_OPEN


template <int dim>
class Quadrature : public Subscriptor
{
public:
  typedef Quadrature<dim-1> SubQuadrature;

  explicit Quadrature (const unsigned int n_quadrature_points = 0);

  Quadrature (const SubQuadrature &,
              const Quadrature<1> &);

  explicit Quadrature (const Quadrature<dim != 1 ? 1 : 0> &quadrature_1d);

  Quadrature (const Quadrature<dim> &q);

  Quadrature (Quadrature<dim> &&) noexcept = default;

  Quadrature (const std::vector<Point<dim> > &points,
              const std::vector<double>      &weights);

  Quadrature (const std::vector<Point<dim> > &points);

  Quadrature (const Point<dim> &point);

  virtual ~Quadrature () = default;

  Quadrature &operator = (const Quadrature<dim> &);

  Quadrature &operator = (Quadrature<dim> &&) = default; // NOLINT

  bool operator == (const Quadrature<dim> &p) const;

  void initialize(const std::vector<Point<dim> > &points,
                  const std::vector<double>      &weights);

  unsigned int size () const;

  const Point<dim> &point (const unsigned int i) const;

  const std::vector<Point<dim> > &get_points () const;

  double weight (const unsigned int i) const;

  const std::vector<double> &get_weights () const;

  std::size_t memory_consumption () const;

  template <class Archive>
  void serialize (Archive &ar, const unsigned int version);

  bool is_tensor_product() const;

  typename std::conditional<dim==1, std::array<Quadrature<1>, dim>,const std::array<Quadrature<1>,dim>&>::type
  get_tensor_basis() const;

protected:
  std::vector<Point<dim> > quadrature_points;

  std::vector<double>      weights;

  bool is_tensor_product_flag;

  std::unique_ptr<std::array<Quadrature<1>,dim>> tensor_basis;
};


template <int dim>
class QAnisotropic : public Quadrature<dim>
{
public:
  QAnisotropic(const Quadrature<1> &qx);

  QAnisotropic(const Quadrature<1> &qx,
               const Quadrature<1> &qy);

  QAnisotropic(const Quadrature<1> &qx,
               const Quadrature<1> &qy,
               const Quadrature<1> &qz);
};


template <int dim>
class QIterated : public Quadrature<dim>
{
public:
  QIterated (const Quadrature<1> &base_quadrature,
             const unsigned int   n_copies);

  DeclExceptionMsg (ExcInvalidQuadratureFormula,
                    "The quadrature formula you provided cannot be used "
                    "as the basis for iteration.");
};



#ifndef DOXYGEN

// -------------------  inline and template functions ----------------


template <int dim>
inline
unsigned int
Quadrature<dim>::size () const
{
  return weights.size();
}


template <int dim>
inline
const Point<dim> &
Quadrature<dim>::point (const unsigned int i) const
{
  AssertIndexRange(i, size());
  return quadrature_points[i];
}



template <int dim>
double
Quadrature<dim>::weight (const unsigned int i) const
{
  AssertIndexRange(i, size());
  return weights[i];
}



template <int dim>
inline
const std::vector<Point<dim> > &
Quadrature<dim>::get_points () const
{
  return quadrature_points;
}



template <int dim>
inline
const std::vector<double> &
Quadrature<dim>::get_weights () const
{
  return weights;
}



template <int dim>
inline
bool
Quadrature<dim>::is_tensor_product () const
{
  return is_tensor_product_flag;
}



template <int dim>
template <class Archive>
inline
void
Quadrature<dim>::serialize (Archive &ar, const unsigned int)
{
  // forward to serialization
  // function in the base class.
  ar   &static_cast<Subscriptor &>(*this);

  ar &quadrature_points &weights;
}



/* -------------- declaration of explicit specializations ------------- */

template <>
Quadrature<0>::Quadrature (const unsigned int);
template <>
Quadrature<0>::Quadrature (const Quadrature<-1> &,
                           const Quadrature<1> &);
template <>
Quadrature<0>::Quadrature (const Quadrature<1> &);

template <>
Quadrature<1>::Quadrature (const Quadrature<0> &,
                           const Quadrature<1> &);

template <>
Quadrature<1>::Quadrature (const Quadrature<0> &);

#endif // DOXYGEN
DEAL_II_NAMESPACE_CLOSE

#endif