quadrature.cc

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2019 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.
//
// ---------------------------------------------------------------------

#include <deal.II/base/geometry_info.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/qprojector.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/std_cxx14/memory.h>
#include <deal.II/base/utilities.h>

#include <cmath>
#include <cstdlib>
#include <iterator>

DEAL_II_NAMESPACE_OPEN


template <>
Quadrature<0>::Quadrature(const unsigned int n_q)
  : quadrature_points(n_q)
  , weights(n_q, 0)
  , is_tensor_product_flag(false)
{}



template <int dim>
Quadrature<dim>::Quadrature(const unsigned int n_q)
  : quadrature_points(n_q, Point<dim>())
  , weights(n_q, 0)
  , is_tensor_product_flag(dim == 1)
{}



template <int dim>
void
Quadrature<dim>::initialize(const std::vector<Point<dim>> &p,
                            const std::vector<double> &    w)
{
  AssertDimension(w.size(), p.size());
  quadrature_points = p;
  weights           = w;
}



template <int dim>
Quadrature<dim>::Quadrature(const std::vector<Point<dim>> &points,
                            const std::vector<double> &    weights)
  : quadrature_points(points)
  , weights(weights)
  , is_tensor_product_flag(dim == 1)
{
  Assert(weights.size() == points.size(),
         ExcDimensionMismatch(weights.size(), points.size()));
}



template <int dim>
Quadrature<dim>::Quadrature(const std::vector<Point<dim>> &points)
  : quadrature_points(points)
  , weights(points.size(), std::numeric_limits<double>::infinity())
  , is_tensor_product_flag(dim == 1)
{
  Assert(weights.size() == points.size(),
         ExcDimensionMismatch(weights.size(), points.size()));
}



template <int dim>
Quadrature<dim>::Quadrature(const Point<dim> &point)
  : quadrature_points(std::vector<Point<dim>>(1, point))
  , weights(std::vector<double>(1, 1.))
  , is_tensor_product_flag(true)
  , tensor_basis(new std::array<Quadrature<1>, dim>())
{
  for (unsigned int i = 0; i < dim; ++i)
    {
      const std::vector<Point<1>> quad_vec_1d(1, Point<1>(point[i]));
      (*tensor_basis)[i] = Quadrature<1>(quad_vec_1d, weights);
    }
}



template <>
Quadrature<1>::Quadrature(const Point<1> &point)
  : quadrature_points(std::vector<Point<1>>(1, point))
  , weights(std::vector<double>(1, 1.))
  , is_tensor_product_flag(true)
{}



template <>
Quadrature<0>::Quadrature(const SubQuadrature &, const Quadrature<1> &)
{
  Assert(false, ExcImpossibleInDim(0));
}



template <int dim>
Quadrature<dim>::Quadrature(const SubQuadrature &q1, const Quadrature<1> &q2)
  : quadrature_points(q1.size() * q2.size())
  , weights(q1.size() * q2.size())
  , is_tensor_product_flag(q1.is_tensor_product())
{
  unsigned int present_index = 0;
  for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
    for (unsigned int i1 = 0; i1 < q1.size(); ++i1)
      {
        // compose coordinates of
        // new quadrature point by tensor
        // product in the last component
        for (unsigned int d = 0; d < dim - 1; ++d)
          quadrature_points[present_index](d) = q1.point(i1)(d);
        quadrature_points[present_index](dim - 1) = q2.point(i2)(0);

        weights[present_index] = q1.weight(i1) * q2.weight(i2);

        ++present_index;
      }

#ifdef DEBUG
  if (size() > 0)
    {
      double sum = 0;
      for (unsigned int i = 0; i < size(); ++i)
        sum += weights[i];
      // we cannot guarantee the sum of weights
      // to be exactly one, but it should be
      // near that.
      Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
    }
#endif

  if (is_tensor_product_flag)
    {
      tensor_basis = std_cxx14::make_unique<std::array<Quadrature<1>, dim>>();
      for (unsigned int i = 0; i < dim - 1; ++i)
        (*tensor_basis)[i] = q1.get_tensor_basis()[i];
      (*tensor_basis)[dim - 1] = q2;
    }
}



template <>
Quadrature<1>::Quadrature(const SubQuadrature &, const Quadrature<1> &q2)
  : quadrature_points(q2.size())
  , weights(q2.size())
  , is_tensor_product_flag(true)
{
  unsigned int present_index = 0;
  for (unsigned int i2 = 0; i2 < q2.size(); ++i2)
    {
      // compose coordinates of
      // new quadrature point by tensor
      // product in the last component
      quadrature_points[present_index](0) = q2.point(i2)(0);

      weights[present_index] = q2.weight(i2);

      ++present_index;
    }

#ifdef DEBUG
  if (size() > 0)
    {
      double sum = 0;
      for (unsigned int i = 0; i < size(); ++i)
        sum += weights[i];
      // we cannot guarantee the sum of weights
      // to be exactly one, but it should be
      // near that.
      Assert((sum > 0.999999) && (sum < 1.000001), ExcInternalError());
    }
#endif
}



template <>
Quadrature<0>::Quadrature(const Quadrature<1> &)
  : Subscriptor()
  ,
  //              quadrature_points(1),
  weights(1, 1.)
  , is_tensor_product_flag(false)
{}


template <>
Quadrature<1>::Quadrature(const Quadrature<0> &)
  : Subscriptor()
{
  // this function should never be
  // called -- this should be the
  // copy constructor in 1d...
  Assert(false, ExcImpossibleInDim(1));
}



template <int dim>
Quadrature<dim>::Quadrature(const Quadrature<dim != 1 ? 1 : 0> &q)
  : Subscriptor()
  , quadrature_points(Utilities::fixed_power<dim>(q.size()))
  , weights(Utilities::fixed_power<dim>(q.size()))
  , is_tensor_product_flag(true)
{
  Assert(dim <= 3, ExcNotImplemented());

  const unsigned int n0 = q.size();
  const unsigned int n1 = (dim > 1) ? n0 : 1;
  const unsigned int n2 = (dim > 2) ? n0 : 1;

  unsigned int k = 0;
  for (unsigned int i2 = 0; i2 < n2; ++i2)
    for (unsigned int i1 = 0; i1 < n1; ++i1)
      for (unsigned int i0 = 0; i0 < n0; ++i0)
        {
          quadrature_points[k](0) = q.point(i0)(0);
          if (dim > 1)
            quadrature_points[k](1) = q.point(i1)(0);
          if (dim > 2)
            quadrature_points[k](2) = q.point(i2)(0);
          weights[k] = q.weight(i0);
          if (dim > 1)
            weights[k] *= q.weight(i1);
          if (dim > 2)
            weights[k] *= q.weight(i2);
          ++k;
        }

  tensor_basis = std_cxx14::make_unique<std::array<Quadrature<1>, dim>>();
  for (unsigned int i = 0; i < dim; ++i)
    (*tensor_basis)[i] = q;
}



template <int dim>
Quadrature<dim>::Quadrature(const Quadrature<dim> &q)
  : Subscriptor()
  , quadrature_points(q.quadrature_points)
  , weights(q.weights)
  , is_tensor_product_flag(q.is_tensor_product_flag)
{
  if (dim > 1 && is_tensor_product_flag)
    tensor_basis =
      std_cxx14::make_unique<std::array<Quadrature<1>, dim>>(*q.tensor_basis);
}



template <int dim>
Quadrature<dim> &
Quadrature<dim>::operator=(const Quadrature<dim> &q)
{
  weights                = q.weights;
  quadrature_points      = q.quadrature_points;
  is_tensor_product_flag = q.is_tensor_product_flag;
  if (dim > 1 && is_tensor_product_flag)
    {
      if (tensor_basis == nullptr)
        tensor_basis = std_cxx14::make_unique<std::array<Quadrature<1>, dim>>(
          *q.tensor_basis);
      else
        *tensor_basis = *q.tensor_basis;
    }
  return *this;
}



template <int dim>
bool
Quadrature<dim>::operator==(const Quadrature<dim> &q) const
{
  return ((quadrature_points == q.quadrature_points) && (weights == q.weights));
}



template <int dim>
std::size_t
Quadrature<dim>::memory_consumption() const
{
  return (MemoryConsumption::memory_consumption(quadrature_points) +
          MemoryConsumption::memory_consumption(weights));
}



template <int dim>
typename std::conditional<dim == 1,
                          std::array<Quadrature<1>, dim>,
                          const std::array<Quadrature<1>, dim> &>::type
Quadrature<dim>::get_tensor_basis() const
{
  Assert(this->is_tensor_product_flag == true,
         ExcMessage("This function only makes sense if "
                    "this object represents a tensor product!"));
  Assert(tensor_basis != nullptr, ExcInternalError());

  return *tensor_basis;
}



template <>
std::array<Quadrature<1>, 1>
Quadrature<1>::get_tensor_basis() const
{
  Assert(this->is_tensor_product_flag == true,
         ExcMessage("This function only makes sense if "
                    "this object represents a tensor product!"));

  return std::array<Quadrature<1>, 1>{{*this}};
}



//---------------------------------------------------------------------------
template <int dim>
QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx)
  : Quadrature<dim>(qx.size())
{
  Assert(dim == 1, ExcImpossibleInDim(dim));
  unsigned int k = 0;
  for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
    {
      this->quadrature_points[k](0) = qx.point(k1)(0);
      this->weights[k++]            = qx.weight(k1);
    }
  Assert(k == this->size(), ExcInternalError());
  this->is_tensor_product_flag = true;
}



template <int dim>
QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx,
                                const Quadrature<1> &qy)
  : Quadrature<dim>(qx.size() * qy.size())
{
  Assert(dim == 2, ExcImpossibleInDim(dim));
}



template <>
QAnisotropic<2>::QAnisotropic(const Quadrature<1> &qx, const Quadrature<1> &qy)
  : Quadrature<2>(qx.size() * qy.size())
{
  unsigned int k = 0;
  for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
    for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
      {
        this->quadrature_points[k](0) = qx.point(k1)(0);
        this->quadrature_points[k](1) = qy.point(k2)(0);
        this->weights[k++]            = qx.weight(k1) * qy.weight(k2);
      }
  Assert(k == this->size(), ExcInternalError());
  this->is_tensor_product_flag = true;
  const std::array<Quadrature<1>, 2> q_array{{qx, qy}};
  this->tensor_basis =
    std_cxx14::make_unique<std::array<Quadrature<1>, 2>>(q_array);
}



template <int dim>
QAnisotropic<dim>::QAnisotropic(const Quadrature<1> &qx,
                                const Quadrature<1> &qy,
                                const Quadrature<1> &qz)
  : Quadrature<dim>(qx.size() * qy.size() * qz.size())
{
  Assert(dim == 3, ExcImpossibleInDim(dim));
}



template <>
QAnisotropic<3>::QAnisotropic(const Quadrature<1> &qx,
                              const Quadrature<1> &qy,
                              const Quadrature<1> &qz)
  : Quadrature<3>(qx.size() * qy.size() * qz.size())
{
  unsigned int k = 0;
  for (unsigned int k3 = 0; k3 < qz.size(); ++k3)
    for (unsigned int k2 = 0; k2 < qy.size(); ++k2)
      for (unsigned int k1 = 0; k1 < qx.size(); ++k1)
        {
          this->quadrature_points[k](0) = qx.point(k1)(0);
          this->quadrature_points[k](1) = qy.point(k2)(0);
          this->quadrature_points[k](2) = qz.point(k3)(0);
          this->weights[k++] = qx.weight(k1) * qy.weight(k2) * qz.weight(k3);
        }
  Assert(k == this->size(), ExcInternalError());
  this->is_tensor_product_flag = true;
  const std::array<Quadrature<1>, 3> q_array{{qx, qy, qz}};
  this->tensor_basis =
    std_cxx14::make_unique<std::array<Quadrature<1>, 3>>(q_array);
}



//---------------------------------------------------------------------------



template <int dim>
Quadrature<2>
QProjector<dim>::reflect(const Quadrature<2> &q)
{
  std::vector<Point<2>> q_points(q.size());
  std::vector<double>   weights(q.size());
  for (unsigned int i = 0; i < q.size(); ++i)
    {
      q_points[i][0] = q.point(i)[1];
      q_points[i][1] = q.point(i)[0];

      weights[i] = q.weight(i);
    }

  return Quadrature<2>(q_points, weights);
}


template <int dim>
Quadrature<2>
QProjector<dim>::rotate(const Quadrature<2> &q, const unsigned int n_times)
{
  std::vector<Point<2>> q_points(q.size());
  std::vector<double>   weights(q.size());
  for (unsigned int i = 0; i < q.size(); ++i)
    {
      switch (n_times % 4)
        {
          case 0:
            // 0 degree
            q_points[i][0] = q.point(i)[0];
            q_points[i][1] = q.point(i)[1];
            break;
          case 1:
            // 90 degree counterclockwise
            q_points[i][0] = 1.0 - q.point(i)[1];
            q_points[i][1] = q.point(i)[0];
            break;
          case 2:
            // 180 degree counterclockwise
            q_points[i][0] = 1.0 - q.point(i)[0];
            q_points[i][1] = 1.0 - q.point(i)[1];
            break;
          case 3:
            // 270 degree counterclockwise
            q_points[i][0] = q.point(i)[1];
            q_points[i][1] = 1.0 - q.point(i)[0];
            break;
        }

      weights[i] = q.weight(i);
    }

  return Quadrature<2>(q_points, weights);
}


template <>
void
QProjector<1>::project_to_face(const Quadrature<0> &,
                               const unsigned int     face_no,
                               std::vector<Point<1>> &q_points)
{
  const unsigned int dim = 1;
  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
  AssertDimension(q_points.size(), 1);

  q_points[0] = Point<dim>(static_cast<double>(face_no));
}



template <>
void
QProjector<2>::project_to_face(const Quadrature<1> &  quadrature,
                               const unsigned int     face_no,
                               std::vector<Point<2>> &q_points)
{
  const unsigned int dim = 2;
  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
  Assert(q_points.size() == quadrature.size(),
         ExcDimensionMismatch(q_points.size(), quadrature.size()));

  for (unsigned int p = 0; p < quadrature.size(); ++p)
    switch (face_no)
      {
        case 0:
          q_points[p] = Point<dim>(0, quadrature.point(p)(0));
          break;
        case 1:
          q_points[p] = Point<dim>(1, quadrature.point(p)(0));
          break;
        case 2:
          q_points[p] = Point<dim>(quadrature.point(p)(0), 0);
          break;
        case 3:
          q_points[p] = Point<dim>(quadrature.point(p)(0), 1);
          break;
        default:
          Assert(false, ExcInternalError());
      }
}



template <>
void
QProjector<3>::project_to_face(const Quadrature<2> &  quadrature,
                               const unsigned int     face_no,
                               std::vector<Point<3>> &q_points)
{
  const unsigned int dim = 3;
  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
  Assert(q_points.size() == quadrature.size(),
         ExcDimensionMismatch(q_points.size(), quadrature.size()));

  for (unsigned int p = 0; p < quadrature.size(); ++p)
    switch (face_no)
      {
        case 0:
          q_points[p] =
            Point<dim>(0, quadrature.point(p)(0), quadrature.point(p)(1));
          break;
        case 1:
          q_points[p] =
            Point<dim>(1, quadrature.point(p)(0), quadrature.point(p)(1));
          break;
        case 2:
          q_points[p] =
            Point<dim>(quadrature.point(p)(1), 0, quadrature.point(p)(0));
          break;
        case 3:
          q_points[p] =
            Point<dim>(quadrature.point(p)(1), 1, quadrature.point(p)(0));
          break;
        case 4:
          q_points[p] =
            Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 0);
          break;
        case 5:
          q_points[p] =
            Point<dim>(quadrature.point(p)(0), quadrature.point(p)(1), 1);
          break;

        default:
          Assert(false, ExcInternalError());
      }
}



template <>
void
QProjector<1>::project_to_subface(const Quadrature<0> &,
                                  const unsigned int face_no,
                                  const unsigned int,
                                  std::vector<Point<1>> &q_points,
                                  const RefinementCase<0> &)
{
  const unsigned int dim = 1;
  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
  AssertDimension(q_points.size(), 1);

  q_points[0] = Point<dim>(static_cast<double>(face_no));
}



template <>
void
QProjector<2>::project_to_subface(const Quadrature<1> &  quadrature,
                                  const unsigned int     face_no,
                                  const unsigned int     subface_no,
                                  std::vector<Point<2>> &q_points,
                                  const RefinementCase<1> &)
{
  const unsigned int dim = 2;
  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
  AssertIndexRange(subface_no, GeometryInfo<dim>::max_children_per_face);

  Assert(q_points.size() == quadrature.size(),
         ExcDimensionMismatch(q_points.size(), quadrature.size()));

  for (unsigned int p = 0; p < quadrature.size(); ++p)
    switch (face_no)
      {
        case 0:
          switch (subface_no)
            {
              case 0:
                q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2);
                break;
              case 1:
                q_points[p] = Point<dim>(0, quadrature.point(p)(0) / 2 + 0.5);
                break;
              default:
                Assert(false, ExcInternalError());
            }
          break;
        case 1:
          switch (subface_no)
            {
              case 0:
                q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2);
                break;
              case 1:
                q_points[p] = Point<dim>(1, quadrature.point(p)(0) / 2 + 0.5);
                break;
              default:
                Assert(false, ExcInternalError());
            }
          break;
        case 2:
          switch (subface_no)
            {
              case 0:
                q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 0);
                break;
              case 1:
                q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 0);
                break;
              default:
                Assert(false, ExcInternalError());
            }
          break;
        case 3:
          switch (subface_no)
            {
              case 0:
                q_points[p] = Point<dim>(quadrature.point(p)(0) / 2, 1);
                break;
              case 1:
                q_points[p] = Point<dim>(quadrature.point(p)(0) / 2 + 0.5, 1);
                break;
              default:
                Assert(false, ExcInternalError());
            }
          break;

        default:
          Assert(false, ExcInternalError());
      }
}



template <>
void
QProjector<3>::project_to_subface(const Quadrature<2> &    quadrature,
                                  const unsigned int       face_no,
                                  const unsigned int       subface_no,
                                  std::vector<Point<3>> &  q_points,
                                  const RefinementCase<2> &ref_case)
{
  const unsigned int dim = 3;
  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
  AssertIndexRange(subface_no, GeometryInfo<dim>::max_children_per_face);
  Assert(q_points.size() == quadrature.size(),
         ExcDimensionMismatch(q_points.size(), quadrature.size()));

  // one coordinate is at a const value. for
  // faces 0, 2 and 4 this value is 0.0, for
  // faces 1, 3 and 5 it is 1.0
  double const_value = face_no % 2;
  // local 2d coordinates are xi and eta,
  // global 3d coordinates are x, y and
  // z. those have to be mapped. the following
  // indices tell, which global coordinate
  // (0->x, 1->y, 2->z) corresponds to which
  // local one
  unsigned int xi_index    = numbers::invalid_unsigned_int,
               eta_index   = numbers::invalid_unsigned_int,
               const_index = face_no / 2;
  // the xi and eta values have to be scaled
  // (by factor 0.5 or factor 1.0) depending on
  // the refinement case and translated (by 0.0
  // or 0.5) depending on the refinement case
  // and subface_no.
  double xi_scale = 1.0, eta_scale = 1.0, xi_translation = 0.0,
         eta_translation = 0.0;
  // set the index mapping between local and
  // global coordinates
  switch (face_no / 2)
    {
      case 0:
        xi_index  = 1;
        eta_index = 2;
        break;
      case 1:
        xi_index  = 2;
        eta_index = 0;
        break;
      case 2:
        xi_index  = 0;
        eta_index = 1;
        break;
    }
  // set the scale and translation parameter
  // for individual subfaces
  switch (ref_case)
    {
      case RefinementCase<dim - 1>::cut_x:
        xi_scale       = 0.5;
        xi_translation = subface_no % 2 * 0.5;
        break;
      case RefinementCase<dim - 1>::cut_y:
        eta_scale       = 0.5;
        eta_translation = subface_no % 2 * 0.5;
        break;
      case RefinementCase<dim - 1>::cut_xy:
        xi_scale        = 0.5;
        eta_scale       = 0.5;
        xi_translation  = int(subface_no % 2) * 0.5;
        eta_translation = int(subface_no / 2) * 0.5;
        break;
      default:
        Assert(false, ExcInternalError());
        break;
    }
  // finally, compute the scaled, translated,
  // projected quadrature points
  for (unsigned int p = 0; p < quadrature.size(); ++p)
    {
      q_points[p][xi_index] =
        xi_scale * quadrature.point(p)(0) + xi_translation;
      q_points[p][eta_index] =
        eta_scale * quadrature.point(p)(1) + eta_translation;
      q_points[p][const_index] = const_value;
    }
}


template <>
Quadrature<1>
QProjector<1>::project_to_all_faces(const Quadrature<0> &quadrature)
{
  const unsigned int dim = 1;

  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell;

  // first fix quadrature points
  std::vector<Point<dim>> q_points;
  q_points.reserve(n_points * n_faces);
  std::vector<Point<dim>> help(n_points);


  // project to each face and append
  // results
  for (unsigned int face = 0; face < n_faces; ++face)
    {
      project_to_face(quadrature, face, help);
      std::copy(help.begin(), help.end(), std::back_inserter(q_points));
    }

  // next copy over weights
  std::vector<double> weights;
  weights.reserve(n_points * n_faces);
  for (unsigned int face = 0; face < n_faces; ++face)
    std::copy(quadrature.get_weights().begin(),
              quadrature.get_weights().end(),
              std::back_inserter(weights));

  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
  Assert(weights.size() == n_points * n_faces, ExcInternalError());

  return Quadrature<dim>(q_points, weights);
}



template <>
Quadrature<2>
QProjector<2>::project_to_all_faces(const SubQuadrature &quadrature)
{
  const unsigned int dim = 2;

  const unsigned int n_points = quadrature.size(),
                     n_faces  = GeometryInfo<dim>::faces_per_cell;

  // first fix quadrature points
  std::vector<Point<dim>> q_points;
  q_points.reserve(n_points * n_faces);
  std::vector<Point<dim>> help(n_points);

  // project to each face and append
  // results
  for (unsigned int face = 0; face < n_faces; ++face)
    {
      project_to_face(quadrature, face, help);
      std::copy(help.begin(), help.end(), std::back_inserter(q_points));
    }

  // next copy over weights
  std::vector<double> weights;
  weights.reserve(n_points * n_faces);
  for (unsigned int face = 0; face < n_faces; ++face)
    std::copy(quadrature.get_weights().begin(),
              quadrature.get_weights().end(),
              std::back_inserter(weights));

  Assert(q_points.size() == n_points * n_faces, ExcInternalError());
  Assert(weights.size() == n_points * n_faces, ExcInternalError());

  return Quadrature<dim>(q_points, weights);
}



template <>
Quadrature<3>
QProjector<3>::project_to_all_faces(const SubQuadrature &quadrature)
{
  const unsigned int dim = 3;

  SubQuadrature q_reflected = reflect(quadrature);
  SubQuadrature q[8]        = {quadrature,
                        rotate(quadrature, 1),
                        rotate(quadrature, 2),
                        rotate(quadrature, 3),
                        q_reflected,
                        rotate(q_reflected, 3),
                        rotate(q_reflected, 2),
                        rotate(q_reflected, 1)};



  const unsigned int n_points = quadrature.size(),
                     n_faces  = GeometryInfo<dim>::faces_per_cell;

  // first fix quadrature points
  std::vector<Point<dim>> q_points;
  q_points.reserve(n_points * n_faces * 8);
  std::vector<Point<dim>> help(n_points);

  std::vector<double> weights;
  weights.reserve(n_points * n_faces * 8);

  // do the following for all possible
  // mutations of a face (mutation==0
  // corresponds to a face with standard
  // orientation, no flip and no rotation)
  for (const auto &mutation : q)
    {
      // project to each face and append
      // results
      for (unsigned int face = 0; face < n_faces; ++face)
        {
          project_to_face(mutation, face, help);
          std::copy(help.begin(), help.end(), std::back_inserter(q_points));
        }

      // next copy over weights
      for (unsigned int face = 0; face < n_faces; ++face)
        std::copy(mutation.get_weights().begin(),
                  mutation.get_weights().end(),
                  std::back_inserter(weights));
    }


  Assert(q_points.size() == n_points * n_faces * 8, ExcInternalError());
  Assert(weights.size() == n_points * n_faces * 8, ExcInternalError());

  return Quadrature<dim>(q_points, weights);
}



template <>
Quadrature<1>
QProjector<1>::project_to_all_subfaces(const Quadrature<0> &quadrature)
{
  const unsigned int dim = 1;

  const unsigned int n_points = 1, n_faces = GeometryInfo<dim>::faces_per_cell,
                     subfaces_per_face =
                       GeometryInfo<dim>::max_children_per_face;

  // first fix quadrature points
  std::vector<Point<dim>> q_points;
  q_points.reserve(n_points * n_faces * subfaces_per_face);
  std::vector<Point<dim>> help(n_points);

  // project to each face and copy
  // results
  for (unsigned int face = 0; face < n_faces; ++face)
    for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
      {
        project_to_subface(quadrature, face, subface, help);
        std::copy(help.begin(), help.end(), std::back_inserter(q_points));
      }

  // next copy over weights
  std::vector<double> weights;
  weights.reserve(n_points * n_faces * subfaces_per_face);
  for (unsigned int face = 0; face < n_faces; ++face)
    for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
      std::copy(quadrature.get_weights().begin(),
                quadrature.get_weights().end(),
                std::back_inserter(weights));

  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
         ExcInternalError());
  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
         ExcInternalError());

  return Quadrature<dim>(q_points, weights);
}



template <>
Quadrature<2>
QProjector<2>::project_to_all_subfaces(const SubQuadrature &quadrature)
{
  const unsigned int dim = 2;

  const unsigned int n_points = quadrature.size(),
                     n_faces  = GeometryInfo<dim>::faces_per_cell,
                     subfaces_per_face =
                       GeometryInfo<dim>::max_children_per_face;

  // first fix quadrature points
  std::vector<Point<dim>> q_points;
  q_points.reserve(n_points * n_faces * subfaces_per_face);
  std::vector<Point<dim>> help(n_points);

  // project to each face and copy
  // results
  for (unsigned int face = 0; face < n_faces; ++face)
    for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
      {
        project_to_subface(quadrature, face, subface, help);
        std::copy(help.begin(), help.end(), std::back_inserter(q_points));
      }

  // next copy over weights
  std::vector<double> weights;
  weights.reserve(n_points * n_faces * subfaces_per_face);
  for (unsigned int face = 0; face < n_faces; ++face)
    for (unsigned int subface = 0; subface < subfaces_per_face; ++subface)
      std::copy(quadrature.get_weights().begin(),
                quadrature.get_weights().end(),
                std::back_inserter(weights));

  Assert(q_points.size() == n_points * n_faces * subfaces_per_face,
         ExcInternalError());
  Assert(weights.size() == n_points * n_faces * subfaces_per_face,
         ExcInternalError());

  return Quadrature<dim>(q_points, weights);
}



template <>
Quadrature<3>
QProjector<3>::project_to_all_subfaces(const SubQuadrature &quadrature)
{
  const unsigned int dim         = 3;
  SubQuadrature      q_reflected = reflect(quadrature);
  SubQuadrature      q[8]        = {quadrature,
                        rotate(quadrature, 1),
                        rotate(quadrature, 2),
                        rotate(quadrature, 3),
                        q_reflected,
                        rotate(q_reflected, 3),
                        rotate(q_reflected, 2),
                        rotate(q_reflected, 1)};

  const unsigned int n_points = quadrature.size(),
                     n_faces  = GeometryInfo<dim>::faces_per_cell,
                     total_subfaces_per_face = 2 + 2 + 4;

  // first fix quadrature points
  std::vector<Point<dim>> q_points;
  q_points.reserve(n_points * n_faces * total_subfaces_per_face * 8);
  std::vector<Point<dim>> help(n_points);

  std::vector<double> weights;
  weights.reserve(n_points * n_faces * total_subfaces_per_face * 8);

  // do the following for all possible
  // mutations of a face (mutation==0
  // corresponds to a face with standard
  // orientation, no flip and no rotation)
  for (const auto &mutation : q)
    {
      // project to each face and copy
      // results
      for (unsigned int face = 0; face < n_faces; ++face)
        for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
             ref_case >= RefinementCase<dim - 1>::cut_x;
             --ref_case)
          for (unsigned int subface = 0;
               subface < GeometryInfo<dim - 1>::n_children(
                           RefinementCase<dim - 1>(ref_case));
               ++subface)
            {
              project_to_subface(mutation,
                                 face,
                                 subface,
                                 help,
                                 RefinementCase<dim - 1>(ref_case));
              std::copy(help.begin(), help.end(), std::back_inserter(q_points));
            }

      // next copy over weights
      for (unsigned int face = 0; face < n_faces; ++face)
        for (unsigned int ref_case = RefinementCase<dim - 1>::cut_xy;
             ref_case >= RefinementCase<dim - 1>::cut_x;
             --ref_case)
          for (unsigned int subface = 0;
               subface < GeometryInfo<dim - 1>::n_children(
                           RefinementCase<dim - 1>(ref_case));
               ++subface)
            std::copy(mutation.get_weights().begin(),
                      mutation.get_weights().end(),
                      std::back_inserter(weights));
    }

  Assert(q_points.size() == n_points * n_faces * total_subfaces_per_face * 8,
         ExcInternalError());
  Assert(weights.size() == n_points * n_faces * total_subfaces_per_face * 8,
         ExcInternalError());

  return Quadrature<dim>(q_points, weights);
}



// This function is not used in the library
template <int dim>
Quadrature<dim>
QProjector<dim>::project_to_child(const Quadrature<dim> &quadrature,
                                  const unsigned int     child_no)
{
  Assert(child_no < GeometryInfo<dim>::max_children_per_cell,
         ExcIndexRange(child_no, 0, GeometryInfo<dim>::max_children_per_cell));

  const unsigned int n_q_points = quadrature.size();

  std::vector<Point<dim>> q_points(n_q_points);
  for (unsigned int i = 0; i < n_q_points; ++i)
    q_points[i] =
      GeometryInfo<dim>::child_to_cell_coordinates(quadrature.point(i),
                                                   child_no);

  // for the weights, things are
  // equally simple: copy them and
  // scale them
  std::vector<double> weights = quadrature.get_weights();
  for (unsigned int i = 0; i < n_q_points; ++i)
    weights[i] *= (1. / GeometryInfo<dim>::max_children_per_cell);

  return Quadrature<dim>(q_points, weights);
}


template <int dim>
Quadrature<dim>
QProjector<dim>::project_to_all_children(const Quadrature<dim> &quadrature)
{
  const unsigned int n_points   = quadrature.size(),
                     n_children = GeometryInfo<dim>::max_children_per_cell;

  std::vector<Point<dim>> q_points(n_points * n_children);
  std::vector<double>     weights(n_points * n_children);

  // project to each child and copy
  // results
  for (unsigned int child = 0; child < n_children; ++child)
    {
      Quadrature<dim> help = project_to_child(quadrature, child);
      for (unsigned int i = 0; i < n_points; ++i)
        {
          q_points[child * n_points + i] = help.point(i);
          weights[child * n_points + i]  = help.weight(i);
        }
    }
  return Quadrature<dim>(q_points, weights);
}



template <int dim>
Quadrature<dim>
QProjector<dim>::project_to_line(const Quadrature<1> &quadrature,
                                 const Point<dim> &   p1,
                                 const Point<dim> &   p2)
{
  const unsigned int      n = quadrature.size();
  std::vector<Point<dim>> points(n);
  std::vector<double>     weights(n);
  const double            length = p1.distance(p2);

  for (unsigned int k = 0; k < n; ++k)
    {
      const double alpha = quadrature.point(k)(0);
      points[k]          = alpha * p2;
      points[k] += (1. - alpha) * p1;
      weights[k] = length * quadrature.weight(k);
    }
  return Quadrature<dim>(points, weights);
}



template <int dim>
typename QProjector<dim>::DataSetDescriptor
QProjector<dim>::DataSetDescriptor::face(const unsigned int face_no,
                                         const bool         face_orientation,
                                         const bool         face_flip,
                                         const bool         face_rotation,
                                         const unsigned int n_quadrature_points)
{
  Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError());

  switch (dim)
    {
      case 1:
      case 2:
        return face_no * n_quadrature_points;


      case 3:
        {
          // in 3d, we have to account for faces that
          // have non-standard face orientation, flip
          // and rotation. thus, we have to store
          // _eight_ data sets per face or subface

          // set up a table with the according offsets
          // for non-standard orientation, first index:
          // face_orientation (standard true=1), second
          // index: face_flip (standard false=0), third
          // index: face_rotation (standard false=0)
          //
          // note, that normally we should use the
          // obvious offsets 0,1,2,3,4,5,6,7. However,
          // prior to the changes enabling flipped and
          // rotated faces, in many places of the
          // library the convention was used, that the
          // first dataset with offset 0 corresponds to
          // a face in standard orientation. therefore
          // we use the offsets 4,5,6,7,0,1,2,3 here to
          // stick to that (implicit) convention
          static const unsigned int offset[2][2][2] = {
            {{4 * GeometryInfo<dim>::faces_per_cell,
              5 * GeometryInfo<dim>::
                    faces_per_cell}, // face_orientation=false; face_flip=false;
                                     // face_rotation=false and true
             {6 * GeometryInfo<dim>::faces_per_cell,
              7 * GeometryInfo<dim>::
                    faces_per_cell}}, // face_orientation=false; face_flip=true;
                                      // face_rotation=false and true
            {{0 * GeometryInfo<dim>::faces_per_cell,
              1 * GeometryInfo<dim>::
                    faces_per_cell}, // face_orientation=true;  face_flip=false;
                                     // face_rotation=false and true
             {2 * GeometryInfo<dim>::faces_per_cell,
              3 * GeometryInfo<dim>::
                    faces_per_cell}}}; // face_orientation=true; face_flip=true;
                                       // face_rotation=false and true

          return (
            (face_no + offset[face_orientation][face_flip][face_rotation]) *
            n_quadrature_points);
        }

      default:
        Assert(false, ExcInternalError());
    }
  return numbers::invalid_unsigned_int;
}



template <>
QProjector<1>::DataSetDescriptor
QProjector<1>::DataSetDescriptor::subface(
  const unsigned int face_no,
  const unsigned int subface_no,
  const bool,
  const bool,
  const bool,
  const unsigned int n_quadrature_points,
  const internal::SubfaceCase<1>)
{
  Assert(face_no < GeometryInfo<1>::faces_per_cell, ExcInternalError());
  Assert(subface_no < GeometryInfo<1>::max_children_per_face,
         ExcInternalError());

  return ((face_no * GeometryInfo<1>::max_children_per_face + subface_no) *
          n_quadrature_points);
}



template <>
QProjector<2>::DataSetDescriptor
QProjector<2>::DataSetDescriptor::subface(
  const unsigned int face_no,
  const unsigned int subface_no,
  const bool,
  const bool,
  const bool,
  const unsigned int n_quadrature_points,
  const internal::SubfaceCase<2>)
{
  Assert(face_no < GeometryInfo<2>::faces_per_cell, ExcInternalError());
  Assert(subface_no < GeometryInfo<2>::max_children_per_face,
         ExcInternalError());

  return ((face_no * GeometryInfo<2>::max_children_per_face + subface_no) *
          n_quadrature_points);
}


template <>
QProjector<3>::DataSetDescriptor
QProjector<3>::DataSetDescriptor::subface(
  const unsigned int             face_no,
  const unsigned int             subface_no,
  const bool                     face_orientation,
  const bool                     face_flip,
  const bool                     face_rotation,
  const unsigned int             n_quadrature_points,
  const internal::SubfaceCase<3> ref_case)
{
  const unsigned int dim = 3;

  Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError());
  Assert(subface_no < GeometryInfo<dim>::max_children_per_face,
         ExcInternalError());

  // As the quadrature points created by
  // QProjector are on subfaces in their
  // "standard location" we have to use a
  // permutation of the equivalent subface
  // number in order to respect face
  // orientation, flip and rotation. The
  // information we need here is exactly the
  // same as the
  // GeometryInfo<3>::child_cell_on_face info
  // for the bottom face (face 4) of a hex, as
  // on this the RefineCase of the cell matches
  // that of the face and the subfaces are
  // numbered in the same way as the child
  // cells.

  // in 3d, we have to account for faces that
  // have non-standard face orientation, flip
  // and rotation. thus, we have to store
  // _eight_ data sets per face or subface
  // already for the isotropic
  // case. Additionally, we have three
  // different refinement cases, resulting in
  // <tt>4 + 2 + 2 = 8</tt> different subfaces
  // for each face.
  const unsigned int total_subfaces_per_face = 8;

  // set up a table with the according offsets
  // for non-standard orientation, first index:
  // face_orientation (standard true=1), second
  // index: face_flip (standard false=0), third
  // index: face_rotation (standard false=0)
  //
  // note, that normally we should use the
  // obvious offsets 0,1,2,3,4,5,6,7. However,
  // prior to the changes enabling flipped and
  // rotated faces, in many places of the
  // library the convention was used, that the
  // first dataset with offset 0 corresponds to
  // a face in standard orientation. therefore
  // we use the offsets 4,5,6,7,0,1,2,3 here to
  // stick to that (implicit) convention
  static const unsigned int orientation_offset[2][2][2] = {
    {// face_orientation=false; face_flip=false; face_rotation=false and true
     {4 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
      5 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
     // face_orientation=false; face_flip=true;  face_rotation=false and true
     {6 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
      7 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}},
    {// face_orientation=true;  face_flip=false; face_rotation=false and true
     {0 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
      1 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face},
     // face_orientation=true;  face_flip=true;  face_rotation=false and true
     {2 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face,
      3 * GeometryInfo<dim>::faces_per_cell * total_subfaces_per_face}}};

  // set up a table with the offsets for a
  // given refinement case respecting the
  // corresponding number of subfaces. the
  // index corresponds to (RefineCase::Type - 1)

  // note, that normally we should use the
  // obvious offsets 0,2,6. However, prior to
  // the implementation of anisotropic
  // refinement, in many places of the library
  // the convention was used, that the first
  // dataset with offset 0 corresponds to a
  // standard (isotropic) face
  // refinement. therefore we use the offsets
  // 6,4,0 here to stick to that (implicit)
  // convention
  static const unsigned int ref_case_offset[3] = {
    6, // cut_x
    4, // cut_y
    0  // cut_xy
  };


  // for each subface of a given FaceRefineCase
  // there is a corresponding equivalent
  // subface number of one of the "standard"
  // RefineCases (cut_x, cut_y, cut_xy). Map
  // the given values to those equivalent
  // ones.

  // first, define an invalid number
  static const unsigned int e = numbers::invalid_unsigned_int;

  static const RefinementCase<dim - 1>
    equivalent_refine_case[internal::SubfaceCase<dim>::case_isotropic + 1]
                          [GeometryInfo<3>::max_children_per_face] = {
                            // case_none. there should be only
                            // invalid values here. However, as
                            // this function is also called (in
                            // tests) for cells which have no
                            // refined faces, use isotropic
                            // refinement instead
                            {RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy},
                            // case_x
                            {RefinementCase<dim - 1>::cut_x,
                             RefinementCase<dim - 1>::cut_x,
                             RefinementCase<dim - 1>::no_refinement,
                             RefinementCase<dim - 1>::no_refinement},
                            // case_x1y
                            {RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_x,
                             RefinementCase<dim - 1>::no_refinement},
                            // case_x2y
                            {RefinementCase<dim - 1>::cut_x,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::no_refinement},
                            // case_x1y2y
                            {RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy},
                            // case_y
                            {RefinementCase<dim - 1>::cut_y,
                             RefinementCase<dim - 1>::cut_y,
                             RefinementCase<dim - 1>::no_refinement,
                             RefinementCase<dim - 1>::no_refinement},
                            // case_y1x
                            {RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_y,
                             RefinementCase<dim - 1>::no_refinement},
                            // case_y2x
                            {RefinementCase<dim - 1>::cut_y,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::no_refinement},
                            // case_y1x2x
                            {RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy},
                            // case_xy (case_isotropic)
                            {RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy,
                             RefinementCase<dim - 1>::cut_xy}};

  static const unsigned int
    equivalent_subface_number[internal::SubfaceCase<dim>::case_isotropic + 1]
                             [GeometryInfo<3>::max_children_per_face] = {
                               // case_none, see above
                               {0, 1, 2, 3},
                               // case_x
                               {0, 1, e, e},
                               // case_x1y
                               {0, 2, 1, e},
                               // case_x2y
                               {0, 1, 3, e},
                               // case_x1y2y
                               {0, 2, 1, 3},
                               // case_y
                               {0, 1, e, e},
                               // case_y1x
                               {0, 1, 1, e},
                               // case_y2x
                               {0, 2, 3, e},
                               // case_y1x2x
                               {0, 1, 2, 3},
                               // case_xy (case_isotropic)
                               {0, 1, 2, 3}};

  // If face-orientation or face_rotation are
  // non-standard, cut_x and cut_y have to be
  // exchanged.
  static const RefinementCase<dim - 1> ref_case_permutation[4] = {
    RefinementCase<dim - 1>::no_refinement,
    RefinementCase<dim - 1>::cut_y,
    RefinementCase<dim - 1>::cut_x,
    RefinementCase<dim - 1>::cut_xy};

  // set a corresponding (equivalent)
  // RefineCase and subface number
  const RefinementCase<dim - 1> equ_ref_case =
    equivalent_refine_case[ref_case][subface_no];
  const unsigned int equ_subface_no =
    equivalent_subface_number[ref_case][subface_no];
  // make sure, that we got a valid subface and RefineCase
  Assert(equ_ref_case != RefinementCase<dim>::no_refinement,
         ExcInternalError());
  Assert(equ_subface_no != e, ExcInternalError());
  // now, finally respect non-standard faces
  const RefinementCase<dim - 1> final_ref_case =
    (face_orientation == face_rotation ? ref_case_permutation[equ_ref_case] :
                                         equ_ref_case);

  // what we have now is the number of
  // the subface in the natural
  // orientation of the *face*. what we
  // need to know is the number of the
  // subface concerning the standard face
  // orientation as seen from the *cell*.

  // this mapping is not trivial, but we
  // have done exactly this stuff in the
  // child_cell_on_face function. in
  // order to reduce the amount of code
  // as well as to make maintaining the
  // functionality easier we want to
  // reuse that information. So we note
  // that on the bottom face (face 4) of
  // a hex cell the local x and y
  // coordinates of the face and the cell
  // coincide, thus also the refinement
  // case of the face corresponds to the
  // refinement case of the cell
  // (ignoring cell refinement along the
  // z direction). Using this knowledge
  // we can (ab)use the
  // child_cell_on_face function to do
  // exactly the transformation we are in
  // need of now
  const unsigned int final_subface_no =
    GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(final_ref_case),
                                          4,
                                          equ_subface_no,
                                          face_orientation,
                                          face_flip,
                                          face_rotation,
                                          equ_ref_case);

  return (((face_no * total_subfaces_per_face +
            ref_case_offset[final_ref_case - 1] + final_subface_no) +
           orientation_offset[face_orientation][face_flip][face_rotation]) *
          n_quadrature_points);
}



template <int dim>
Quadrature<dim>
QProjector<dim>::project_to_face(const SubQuadrature &quadrature,
                                 const unsigned int   face_no)
{
  std::vector<Point<dim>> points(quadrature.size());
  project_to_face(quadrature, face_no, points);
  return Quadrature<dim>(points, quadrature.get_weights());
}


template <int dim>
Quadrature<dim>
QProjector<dim>::project_to_subface(const SubQuadrature &          quadrature,
                                    const unsigned int             face_no,
                                    const unsigned int             subface_no,
                                    const RefinementCase<dim - 1> &ref_case)
{
  std::vector<Point<dim>> points(quadrature.size());
  project_to_subface(quadrature, face_no, subface_no, points, ref_case);
  return Quadrature<dim>(points, quadrature.get_weights());
}


// ------------------------------------------------------------ //

namespace internal
{
  namespace QIteratedImplementation
  {
    namespace
    {
      bool
      uses_both_endpoints(const Quadrature<1> &base_quadrature)
      {
        bool at_left = false, at_right = false;
        for (unsigned int i = 0; i < base_quadrature.size(); ++i)
          {
            if (base_quadrature.point(i) == Point<1>(0.0))
              at_left = true;
            if (base_quadrature.point(i) == Point<1>(1.0))
              at_right = true;
          }

        return (at_left && at_right);
      }
    } // namespace
  }   // namespace QIteratedImplementation
} // namespace internal

// template <>
// void
// QIterated<1>::fill(Quadrature<1>& dst,
//                 const Quadrature<1> &base_quadrature,
//                 const unsigned int   n_copies)
// {
//   Assert (n_copies > 0, ExcZero());
//   Assert (base_quadrature.size() > 0, ExcZero());

//   const unsigned int np =
//     uses_both_endpoints(base_quadrature)
//     ? (base_quadrature.size()-1) * n_copies + 1
//     : base_quadrature.size() * n_copies;

//   dst.quadrature_points.resize(np);
//   dst.weights.resize(np);

//   if (!uses_both_endpoints(base_quadrature))
//                                   // we don't have to skip some
//                                   // points in order to get a
//                                   // reasonable quadrature formula
//     {
//       unsigned int next_point = 0;
//       for (unsigned int copy=0; copy<n_copies; ++copy)
//      for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
//        {
//          dst.quadrature_points[next_point](0)
//            = (copy + base_quadrature.point(q_point)(0)) / n_copies;
//          dst.weights[next_point]
//            = base_quadrature.weight(q_point) / n_copies;
//          ++next_point;
//        }
//     }
//   else
//                                   // skip doubly available points
//     {
//       unsigned int next_point = 0;

//                                     // first find out the weights of
//                                     // the left and the right boundary
//                                     // points. note that these usually
//                                     // are but need not necessarily be
//                                     // the same
//       double double_point_weight = 0;
//       unsigned int n_end_points = 0;
//       for (unsigned int i=0; i<base_quadrature.size(); ++i)
//                                       // add up the weight if this
//                                       // is an endpoint
//      if ((base_quadrature.point(i)(0) == 0.) ||
//          (base_quadrature.point(i)(0) == 1.))
//        {
//          double_point_weight += base_quadrature.weight(i);
//          ++n_end_points;
//        }
//                                     // scale the weight correctly
//       double_point_weight /= n_copies;

//                                     // make sure the base quadrature formula
//                                     // has only one quadrature point
//                                     // per end point
//       Assert (n_end_points == 2, ExcInvalidQuadratureFormula());


//       for (unsigned int copy=0; copy<n_copies; ++copy)
//      for (unsigned int q_point=0; q_point<base_quadrature.size(); ++q_point)
//        {
//                                           // skip the left point of
//                                           // this copy since we
//                                           // have already entered
//                                           // it the last time
//          if ((copy > 0) &&
//              (base_quadrature.point(q_point)(0) == 0.))
//            continue;

//          dst.quadrature_points[next_point](0)
//            = (copy+base_quadrature.point(q_point)(0)) / n_copies;

//                                           // if this is the
//                                           // rightmost point of one
//                                           // of the non-last
//                                           // copies: give it the
//                                           // double weight
//          if ((copy != n_copies-1) &&
//              (base_quadrature.point(q_point)(0) == 1.))
//            dst.weights[next_point] = double_point_weight;
//          else
//            dst.weights[next_point] = base_quadrature.weight(q_point) /
//                                        n_copies;

//          ++next_point;
//        }
//     }

// #if DEBUG
//   double sum_of_weights = 0;
//   for (unsigned int i=0; i<dst.size(); ++i)
//     sum_of_weights += dst.weight(i);
//   Assert (std::fabs(sum_of_weights-1) < 1e-15,
//        ExcInternalError());
// #endif

// }


template <>
QIterated<0>::QIterated(const Quadrature<1> &, const unsigned int)
  : Quadrature<0>()
{
  Assert(false, ExcNotImplemented());
}



template <>
QIterated<1>::QIterated(const Quadrature<1> &base_quadrature,
                        const unsigned int   n_copies)
  : Quadrature<1>(
      internal::QIteratedImplementation::uses_both_endpoints(base_quadrature) ?
        (base_quadrature.size() - 1) * n_copies + 1 :
        base_quadrature.size() * n_copies)
{
  //  fill(*this, base_quadrature, n_copies);
  Assert(base_quadrature.size() > 0, ExcNotInitialized());
  Assert(n_copies > 0, ExcZero());

  if (!internal::QIteratedImplementation::uses_both_endpoints(base_quadrature))
    // we don't have to skip some
    // points in order to get a
    // reasonable quadrature formula
    {
      unsigned int next_point = 0;
      for (unsigned int copy = 0; copy < n_copies; ++copy)
        for (unsigned int q_point = 0; q_point < base_quadrature.size();
             ++q_point)
          {
            this->quadrature_points[next_point] =
              Point<1>(base_quadrature.point(q_point)(0) / n_copies +
                       (1.0 * copy) / n_copies);
            this->weights[next_point] =
              base_quadrature.weight(q_point) / n_copies;

            ++next_point;
          }
    }
  else
    // skip doubly available points
    {
      unsigned int next_point = 0;

      // first find out the weights of
      // the left and the right boundary
      // points. note that these usually
      // are but need not necessarily be
      // the same
      double       double_point_weight = 0;
      unsigned int n_end_points        = 0;
      for (unsigned int i = 0; i < base_quadrature.size(); ++i)
        // add up the weight if this
        // is an endpoint
        if ((base_quadrature.point(i) == Point<1>(0.0)) ||
            (base_quadrature.point(i) == Point<1>(1.0)))
          {
            double_point_weight += base_quadrature.weight(i);
            ++n_end_points;
          }
      // scale the weight correctly
      double_point_weight /= n_copies;

      // make sure the base quadrature formula
      // has only one quadrature point
      // per end point
      Assert(n_end_points == 2, ExcInvalidQuadratureFormula());


      for (unsigned int copy = 0; copy < n_copies; ++copy)
        for (unsigned int q_point = 0; q_point < base_quadrature.size();
             ++q_point)
          {
            // skip the left point of
            // this copy since we
            // have already entered
            // it the last time
            if ((copy > 0) && (base_quadrature.point(q_point) == Point<1>(0.0)))
              continue;

            this->quadrature_points[next_point] =
              Point<1>(base_quadrature.point(q_point)(0) / n_copies +
                       (1.0 * copy) / n_copies);

            // if this is the
            // rightmost point of one
            // of the non-last
            // copies: give it the
            // double weight
            if ((copy != n_copies - 1) &&
                (base_quadrature.point(q_point) == Point<1>(1.0)))
              this->weights[next_point] = double_point_weight;
            else
              this->weights[next_point] =
                base_quadrature.weight(q_point) / n_copies;

            ++next_point;
          }
    }

#if DEBUG
  double sum_of_weights = 0;
  for (unsigned int i = 0; i < this->size(); ++i)
    sum_of_weights += this->weight(i);
  Assert(std::fabs(sum_of_weights - 1) < 1e-13, ExcInternalError());
#endif
}


// template <int dim>
// void
// QIterated<dim>::fill(Quadrature<dim>&, const Quadrature<1>&, unsigned int)
// {
//   Assert(false, ExcNotImplemented());
// }


// construct higher dimensional quadrature formula by tensor product
// of lower dimensional iterated quadrature formulae
template <int dim>
QIterated<dim>::QIterated(const Quadrature<1> &base_quadrature,
                          const unsigned int   N)
  : Quadrature<dim>(QIterated<dim - 1>(base_quadrature, N),
                    QIterated<1>(base_quadrature, N))
{}



// explicit instantiations; note: we need them all for all dimensions
template class Quadrature<0>;
template class Quadrature<1>;
template class Quadrature<2>;
template class Quadrature<3>;
template class QAnisotropic<1>;
template class QAnisotropic<2>;
template class QAnisotropic<3>;
template class QIterated<1>;
template class QIterated<2>;
template class QIterated<3>;
template class QProjector<1>;
template class QProjector<2>;
template class QProjector<3>;

DEAL_II_NAMESPACE_CLOSE