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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#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 <array>
#include <memory>
#include <vector>

DEAL_II_NAMESPACE_OPEN


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

  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() override = 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