#include <deal.II/base/geometry_info.h>
#include <deal.II/base/quadrature.h>
+#include <deal.II/grid/reference_cell.h>
+
DEAL_II_NAMESPACE_OPEN
* Compute the quadrature points on the cell if the given quadrature formula
* is used on face <tt>face_no</tt>. For further details, see the general
* doc for this class.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static void
+ DEAL_II_DEPRECATED static void
project_to_face(const SubQuadrature & quadrature,
const unsigned int face_no,
std::vector<Point<dim>> &q_points);
+ /**
+ * Compute the quadrature points on the cell if the given quadrature formula
+ * is used on face <tt>face_no</tt>. For further details, see the general
+ * doc for this class.
+ */
+ static void
+ project_to_face(const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature,
+ const unsigned int face_no,
+ std::vector<Point<dim>> & q_points);
+
/**
* Compute the cell quadrature formula corresponding to using
* <tt>quadrature</tt> on face <tt>face_no</tt>. For further details, see
* the general doc for this class.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static Quadrature<dim>
+ DEAL_II_DEPRECATED static Quadrature<dim>
project_to_face(const SubQuadrature &quadrature, const unsigned int face_no);
+ /**
+ * Compute the cell quadrature formula corresponding to using
+ * <tt>quadrature</tt> on face <tt>face_no</tt>. For further details, see
+ * the general doc for this class.
+ */
+ static Quadrature<dim>
+ project_to_face(const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature,
+ const unsigned int face_no);
+
/**
* Compute the quadrature points on the cell if the given quadrature formula
* is used on face <tt>face_no</tt>, subface number <tt>subface_no</tt>
*
* @note Only the points are transformed. The quadrature weights are the
* same as those of the original rule.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static void
+ DEAL_II_DEPRECATED static void
project_to_subface(const SubQuadrature & quadrature,
const unsigned int face_no,
const unsigned int subface_no,
const RefinementCase<dim - 1> &ref_case =
RefinementCase<dim - 1>::isotropic_refinement);
+ /**
+ * Compute the quadrature points on the cell if the given quadrature formula
+ * is used on face <tt>face_no</tt>, subface number <tt>subface_no</tt>
+ * corresponding to RefineCase::Type <tt>ref_case</tt>. The last argument is
+ * only used in 3D.
+ *
+ * @note Only the points are transformed. The quadrature weights are the
+ * same as those of the original rule.
+ */
+ static void
+ project_to_subface(const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ std::vector<Point<dim>> & q_points,
+ const RefinementCase<dim - 1> &ref_case =
+ RefinementCase<dim - 1>::isotropic_refinement);
+
/**
* Compute the cell quadrature formula corresponding to using
* <tt>quadrature</tt> on subface <tt>subface_no</tt> of face
*
* @note Only the points are transformed. The quadrature weights are the
* same as those of the original rule.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static Quadrature<dim>
+ DEAL_II_DEPRECATED static Quadrature<dim>
project_to_subface(const SubQuadrature & quadrature,
const unsigned int face_no,
const unsigned int subface_no,
const RefinementCase<dim - 1> &ref_case =
RefinementCase<dim - 1>::isotropic_refinement);
+ /**
+ * Compute the cell quadrature formula corresponding to using
+ * <tt>quadrature</tt> on subface <tt>subface_no</tt> of face
+ * <tt>face_no</tt> with RefinementCase<dim-1> <tt>ref_case</tt>. The last
+ * argument is only used in 3D.
+ *
+ * @note Only the points are transformed. The quadrature weights are the
+ * same as those of the original rule.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
+ */
+ static Quadrature<dim>
+ project_to_subface(const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const RefinementCase<dim - 1> &ref_case =
+ RefinementCase<dim - 1>::isotropic_refinement);
+
/**
* Take a face quadrature formula and generate a cell quadrature formula
* from it where the quadrature points of the given argument are projected
*
* @note In 3D, this function produces eight sets of quadrature points for
* each face, in order to cope possibly different orientations of the mesh.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static Quadrature<dim>
+ DEAL_II_DEPRECATED static Quadrature<dim>
project_to_all_faces(const SubQuadrature &quadrature);
+ /**
+ * Take a face quadrature formula and generate a cell quadrature formula
+ * from it where the quadrature points of the given argument are projected
+ * on all faces.
+ *
+ * The weights of the new rule are replications of the original weights.
+ * Thus, the sum of the weights is not one, but the number of faces, which
+ * is the surface of the reference cell.
+ *
+ * This in particular allows us to extract a subset of points corresponding
+ * to a single face and use it as a quadrature on this face, as is done in
+ * FEFaceValues.
+ *
+ * @note In 3D, this function produces eight sets of quadrature points for
+ * each face, in order to cope possibly different orientations of the mesh.
+ */
+ static Quadrature<dim>
+ project_to_all_faces(const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature);
+
/**
* Take a face quadrature formula and generate a cell quadrature formula
* from it where the quadrature points of the given argument are projected
* This in particular allows us to extract a subset of points corresponding
* to a single subface and use it as a quadrature on this face, as is done
* in FESubfaceValues.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static Quadrature<dim>
+ DEAL_II_DEPRECATED static Quadrature<dim>
project_to_all_subfaces(const SubQuadrature &quadrature);
+ /**
+ * Take a face quadrature formula and generate a cell quadrature formula
+ * from it where the quadrature points of the given argument are projected
+ * on all subfaces.
+ *
+ * Like in project_to_all_faces(), the weights of the new rule sum up to the
+ * number of faces (not subfaces), which is the surface of the reference
+ * cell.
+ *
+ * This in particular allows us to extract a subset of points corresponding
+ * to a single subface and use it as a quadrature on this face, as is done
+ * in FESubfaceValues.
+ */
+ static Quadrature<dim>
+ project_to_all_subfaces(const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature);
+
/**
* Project a given quadrature formula to a child of a cell. You may want to
* use this function in case you want to extend an integral only over the
* As integration using this quadrature formula now only extends over a
* fraction of the cell, the weights of the resulting object are divided by
* GeometryInfo<dim>::children_per_cell.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static Quadrature<dim>
+ DEAL_II_DEPRECATED static Quadrature<dim>
project_to_child(const Quadrature<dim> &quadrature,
const unsigned int child_no);
+ /**
+ * Project a given quadrature formula to a child of a cell. You may want to
+ * use this function in case you want to extend an integral only over the
+ * area which a potential child would occupy. The child numbering is the
+ * same as the children would be numbered upon refinement of the cell.
+ *
+ * As integration using this quadrature formula now only extends over a
+ * fraction of the cell, the weights of the resulting object are divided by
+ * GeometryInfo<dim>::children_per_cell.
+ */
+ static Quadrature<dim>
+ project_to_child(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<dim> & quadrature,
+ const unsigned int child_no);
+
/**
* Project a quadrature rule to all children of a cell. Similarly to
* project_to_all_subfaces(), this function replicates the formula generated
*
* The child numbering is the same as the children would be numbered upon
* refinement of the cell.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static Quadrature<dim>
+ DEAL_II_DEPRECATED static Quadrature<dim>
project_to_all_children(const Quadrature<dim> &quadrature);
+ /**
+ * Project a quadrature rule to all children of a cell. Similarly to
+ * project_to_all_subfaces(), this function replicates the formula generated
+ * by project_to_child() for all children, such that the weights sum up to
+ * one, the volume of the total cell again.
+ *
+ * The child numbering is the same as the children would be numbered upon
+ * refinement of the cell.
+ */
+ static Quadrature<dim>
+ project_to_all_children(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<dim> & quadrature);
+
/**
* Project the one dimensional rule <tt>quadrature</tt> to the straight line
* connecting the points <tt>p1</tt> and <tt>p2</tt>.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static Quadrature<dim>
+ DEAL_II_DEPRECATED static Quadrature<dim>
project_to_line(const Quadrature<1> &quadrature,
const Point<dim> & p1,
const Point<dim> & p2);
+ /**
+ * Project the one dimensional rule <tt>quadrature</tt> to the straight line
+ * connecting the points <tt>p1</tt> and <tt>p2</tt>.
+ */
+ static Quadrature<dim>
+ project_to_line(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<1> & quadrature,
+ const Point<dim> & p1,
+ const Point<dim> & p2);
+
/**
* Since the project_to_all_faces() and project_to_all_subfaces() functions
* chain together the quadrature points and weights of all projections of a
* The last argument denotes the number of quadrature points the lower-
* dimensional face quadrature formula (the one that has been projected
* onto the faces) has.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static DataSetDescriptor
+ DEAL_II_DEPRECATED static 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);
+ /**
+ * Static function to generate an offset object for a given face of a cell
+ * with the given face orientation, flip and rotation. This function of
+ * course is only allowed if <tt>dim>=2</tt>, and the face orientation,
+ * flip and rotation are ignored if the space dimension equals 2.
+ *
+ * The last argument denotes the number of quadrature points the lower-
+ * dimensional face quadrature formula (the one that has been projected
+ * onto the faces) has.
+ */
+ static DataSetDescriptor
+ face(const ReferenceCell::Type reference_cell_type,
+ const unsigned int face_no,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation,
+ const unsigned int n_quadrature_points);
+
/**
* Static function to generate an offset object for a given subface of a
* cell with the given face orientation, flip and rotation. This function
* projected onto the faces) has.
*
* Through the last argument anisotropic refinement can be respected.
+ *
+ * @deprecated This function makes an implicit assumption that the cell is
+ * a line (1D), a quad (2D), or a hex (3D). Use the other version of
+ * this function that takes the reference cell type instead.
*/
- static DataSetDescriptor
+ DEAL_II_DEPRECATED static DataSetDescriptor
subface(const unsigned int face_no,
const unsigned int subface_no,
const bool face_orientation,
const internal::SubfaceCase<dim> ref_case =
internal::SubfaceCase<dim>::case_isotropic);
+ /**
+ * Static function to generate an offset object for a given subface of a
+ * cell with the given face orientation, flip and rotation. This function
+ * of course is only allowed if <tt>dim>=2</tt>, and the face orientation,
+ * flip and rotation are ignored if the space dimension equals 2.
+ *
+ * The last but one argument denotes the number of quadrature points the
+ * lower-dimensional face quadrature formula (the one that has been
+ * projected onto the faces) has.
+ *
+ * Through the last argument anisotropic refinement can be respected.
+ */
+ static DataSetDescriptor
+ subface(const ReferenceCell::Type reference_cell_type,
+ 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<dim> ref_case =
+ internal::SubfaceCase<dim>::case_isotropic);
+
/**
* Conversion operator to an integer denoting the offset of the first
* element of this dataset in the set of quadrature formulas all projected
std::vector<Point<1>> &);
template <>
void
+QProjector<1>::project_to_face(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<0> &,
+ const unsigned int,
+ std::vector<Point<1>> &);
+template <>
+void
QProjector<2>::project_to_face(const Quadrature<1> & quadrature,
const unsigned int face_no,
std::vector<Point<2>> &q_points);
template <>
void
+QProjector<2>::project_to_face(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<1> & quadrature,
+ const unsigned int face_no,
+ std::vector<Point<2>> & q_points);
+template <>
+void
QProjector<3>::project_to_face(const Quadrature<2> & quadrature,
const unsigned int face_no,
std::vector<Point<3>> &q_points);
+template <>
+void
+QProjector<3>::project_to_face(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<2> & quadrature,
+ const unsigned int face_no,
+ std::vector<Point<3>> & q_points);
template <>
Quadrature<1>
QProjector<1>::project_to_all_faces(const Quadrature<0> &quadrature);
+template <>
+Quadrature<1>
+QProjector<1>::project_to_all_faces(
+ const ReferenceCell::Type reference_cell_type,
+ const Quadrature<0> & quadrature);
template <>
const RefinementCase<0> &);
template <>
void
+QProjector<1>::project_to_subface(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<0> &,
+ const unsigned int,
+ const unsigned int,
+ std::vector<Point<1>> &,
+ const RefinementCase<0> &);
+template <>
+void
QProjector<2>::project_to_subface(const Quadrature<1> & quadrature,
const unsigned int face_no,
const unsigned int subface_no,
const RefinementCase<1> &);
template <>
void
+QProjector<2>::project_to_subface(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<1> & quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ std::vector<Point<2>> & q_points,
+ const RefinementCase<1> &);
+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> &face_ref_case);
+template <>
+void
+QProjector<3>::project_to_subface(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<2> & quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ std::vector<Point<3>> & q_points,
+ const RefinementCase<2> & face_ref_case);
template <>
Quadrature<1>
QProjector<1>::project_to_all_subfaces(const Quadrature<0> &quadrature);
+template <>
+Quadrature<1>
+QProjector<1>::project_to_all_subfaces(
+ const ReferenceCell::Type reference_cell_type,
+ const Quadrature<0> & quadrature);
#endif // DOXYGEN
// to take (all data sets for all
// faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset =
- QProjector<dim>::DataSetDescriptor::face(face_no,
+ const auto offset =
+ QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
+ face_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
cell->face_rotation(face_no),
// to take (all data sets for all
// sub-faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset =
- QProjector<dim>::DataSetDescriptor::subface(face_no,
+ const auto offset =
+ QProjector<dim>::DataSetDescriptor::subface(this->reference_cell_type(),
+ face_no,
sub_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
return get_face_data(update_flags,
mapping,
QProjector<dim - 1>::project_to_all_children(
- quadrature),
+ ReferenceCell::get_hypercube(dim - 1), quadrature),
output_data);
}
// hating the anisotropic implementation
QGauss<dim - 1> q_gauss(degree + 1);
const Quadrature<dim> q_fine =
- QProjector<dim>::project_to_face(q_gauss, face_fine);
+ QProjector<dim>::project_to_face(fe.reference_cell_type(),
+ q_gauss,
+ face_fine);
const unsigned int nq = q_fine.size();
FEValues<dim> fine(mapping,
cell_number < GeometryInfo<dim>::max_children_per_face;
++cell_number)
{
- const Quadrature<dim> q_coarse =
- QProjector<dim>::project_to_subface(q_gauss,
- face_coarse,
- cell_number);
+ const Quadrature<dim> q_coarse = QProjector<dim>::project_to_subface(
+ fe.reference_cell_type(), q_gauss, face_coarse, cell_number);
FEValues<dim> coarse(mapping, fe, q_coarse, update_values);
typename Triangulation<dim, spacedim>::active_cell_iterator fine_cell =
Invalid = static_cast<std::uint8_t>(-1)
};
+ /**
+ * Return the correct simplex reference cell type for the given dimension
+ * @p dim.
+ */
+ inline Type
+ get_simplex(const unsigned int dim)
+ {
+ switch (dim)
+ {
+ case 0:
+ return Type::Vertex;
+ case 1:
+ return Type::Line;
+ case 2:
+ return Type::Tri;
+ case 3:
+ return Type::Tet;
+ default:
+ Assert(false, ExcNotImplemented());
+ return Type::Invalid;
+ }
+ }
+
+ /**
+ * Return the correct hypercube reference cell type for the given dimension
+ * @p dim.
+ */
+ inline Type
+ get_hypercube(const unsigned int dim)
+ {
+ switch (dim)
+ {
+ case 0:
+ return Type::Vertex;
+ case 1:
+ return Type::Line;
+ case 2:
+ return Type::Quad;
+ case 3:
+ return Type::Hex;
+ default:
+ Assert(false, ExcNotImplemented());
+ return Type::Invalid;
+ }
+ }
+
namespace internal
{
/**
polynomials_raviart_thomas.cc
polynomials_rt_bubbles.cc
process_grid.cc
+ qprojector.cc
quadrature.cc
quadrature_lib.cc
quadrature_selector.cc
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 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/qprojector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+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> & quadrature,
+ const unsigned int face_no,
+ std::vector<Point<1>> &q_points)
+{
+ project_to_face(ReferenceCell::Type::Line, quadrature, face_no, q_points);
+}
+
+
+
+template <>
+void
+QProjector<1>::project_to_face(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<0> &,
+ const unsigned int face_no,
+ std::vector<Point<1>> &q_points)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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)
+{
+ project_to_face(ReferenceCell::Type::Quad, quadrature, face_no, q_points);
+}
+
+
+
+template <>
+void
+QProjector<2>::project_to_face(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<1> & quadrature,
+ const unsigned int face_no,
+ std::vector<Point<2>> & q_points)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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)
+{
+ project_to_face(ReferenceCell::Type::Hex, quadrature, face_no, q_points);
+}
+
+
+
+template <>
+void
+QProjector<3>::project_to_face(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<2> & quadrature,
+ const unsigned int face_no,
+ std::vector<Point<3>> & q_points)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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> & quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ std::vector<Point<1>> & q_points,
+ const RefinementCase<0> &ref_case)
+{
+ project_to_subface(ReferenceCell::Type::Line,
+ quadrature,
+ face_no,
+ subface_no,
+ q_points,
+ ref_case);
+}
+
+
+
+template <>
+void
+QProjector<1>::project_to_subface(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<0> &,
+ const unsigned int face_no,
+ const unsigned int,
+ std::vector<Point<1>> &q_points,
+ const RefinementCase<0> &)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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> &ref_case)
+{
+ project_to_subface(ReferenceCell::Type::Quad,
+ quadrature,
+ face_no,
+ subface_no,
+ q_points,
+ ref_case);
+}
+
+
+
+template <>
+void
+QProjector<2>::project_to_subface(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<1> & quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ std::vector<Point<2>> & q_points,
+ const RefinementCase<1> &)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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)
+{
+ project_to_subface(ReferenceCell::Type::Hex,
+ quadrature,
+ face_no,
+ subface_no,
+ q_points,
+ ref_case);
+}
+
+
+
+template <>
+void
+QProjector<3>::project_to_subface(const ReferenceCell::Type reference_cell_type,
+ 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)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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)
+{
+ return project_to_all_faces(ReferenceCell::Type::Line, quadrature);
+}
+
+
+template <>
+Quadrature<1>
+QProjector<1>::project_to_all_faces(
+ const ReferenceCell::Type reference_cell_type,
+ const Quadrature<0> & quadrature)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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 ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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<2>
+QProjector<2>::project_to_all_faces(const SubQuadrature &quadrature)
+{
+ return project_to_all_faces(ReferenceCell::Type::Quad, quadrature);
+}
+
+
+
+template <>
+Quadrature<3>
+QProjector<3>::project_to_all_faces(
+ const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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<3>
+QProjector<3>::project_to_all_faces(const SubQuadrature &quadrature)
+{
+ return project_to_all_faces(ReferenceCell::Type::Hex, quadrature);
+}
+
+
+
+template <>
+Quadrature<1>
+QProjector<1>::project_to_all_subfaces(const Quadrature<0> &quadrature)
+{
+ return project_to_all_subfaces(ReferenceCell::Type::Line, quadrature);
+}
+
+
+
+template <>
+Quadrature<1>
+QProjector<1>::project_to_all_subfaces(
+ const ReferenceCell::Type reference_cell_type,
+ const Quadrature<0> & quadrature)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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 ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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<2>
+QProjector<2>::project_to_all_subfaces(const SubQuadrature &quadrature)
+{
+ return project_to_all_subfaces(ReferenceCell::Type::Quad, quadrature);
+}
+
+
+
+template <>
+Quadrature<3>
+QProjector<3>::project_to_all_subfaces(
+ const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature)
+{
+ Assert(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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);
+}
+
+
+
+template <>
+Quadrature<3>
+QProjector<3>::project_to_all_subfaces(const SubQuadrature &quadrature)
+{
+ return project_to_all_subfaces(ReferenceCell::Type::Hex, quadrature);
+}
+
+
+
+// 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)
+{
+ return project_to_child(ReferenceCell::get_hypercube(dim),
+ quadrature,
+ child_no);
+}
+
+
+
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_child(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<dim> & quadrature,
+ const unsigned int child_no)
+{
+ Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
+ ExcNotImplemented());
+ (void)reference_cell_type;
+
+ AssertIndexRange(child_no, 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)
+{
+ return project_to_all_children(ReferenceCell::get_hypercube(dim), quadrature);
+}
+
+
+
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_all_children(
+ const ReferenceCell::Type reference_cell_type,
+ const Quadrature<dim> & quadrature)
+{
+ Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
+ ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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)
+{
+ return project_to_line(ReferenceCell::get_hypercube(dim), quadrature, p1, p2);
+}
+
+
+
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_line(const ReferenceCell::Type reference_cell_type,
+ const Quadrature<1> & quadrature,
+ const Point<dim> & p1,
+ const Point<dim> & p2)
+{
+ Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
+ ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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)
+{
+ return face(ReferenceCell::get_hypercube(dim),
+ face_no,
+ face_orientation,
+ face_flip,
+ face_rotation,
+ n_quadrature_points);
+}
+
+
+
+template <int dim>
+typename QProjector<dim>::DataSetDescriptor
+QProjector<dim>::DataSetDescriptor::face(
+ const ReferenceCell::Type reference_cell_type,
+ const unsigned int face_no,
+ const bool face_orientation,
+ const bool face_flip,
+ const bool face_rotation,
+ const unsigned int n_quadrature_points)
+{
+ Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
+ ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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 ReferenceCell::Type reference_cell_type,
+ 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(reference_cell_type == ReferenceCell::Type::Line, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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<1>::DataSetDescriptor
+QProjector<1>::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<1> ref_case)
+{
+ return subface(ReferenceCell::Type::Line,
+ face_no,
+ subface_no,
+ face_orientation,
+ face_flip,
+ face_rotation,
+ n_quadrature_points,
+ ref_case);
+}
+
+
+
+template <>
+QProjector<2>::DataSetDescriptor
+QProjector<2>::DataSetDescriptor::subface(
+ const ReferenceCell::Type reference_cell_type,
+ 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(reference_cell_type == ReferenceCell::Type::Quad, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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<2>::DataSetDescriptor
+QProjector<2>::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<2> ref_case)
+{
+ return subface(ReferenceCell::Type::Quad,
+ face_no,
+ subface_no,
+ face_orientation,
+ face_flip,
+ face_rotation,
+ n_quadrature_points,
+ ref_case);
+}
+
+
+template <>
+QProjector<3>::DataSetDescriptor
+QProjector<3>::DataSetDescriptor::subface(
+ const ReferenceCell::Type reference_cell_type,
+ 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(reference_cell_type == ReferenceCell::Type::Hex, ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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 <>
+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)
+{
+ return subface(ReferenceCell::Type::Hex,
+ face_no,
+ subface_no,
+ face_orientation,
+ face_flip,
+ face_rotation,
+ n_quadrature_points,
+ ref_case);
+}
+
+
+
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_face(const SubQuadrature &quadrature,
+ const unsigned int face_no)
+{
+ return project_to_face(ReferenceCell::get_hypercube(dim),
+ quadrature,
+ face_no);
+}
+
+
+
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_face(const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature,
+ const unsigned int face_no)
+{
+ Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
+ ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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)
+{
+ return project_to_subface(ReferenceCell::get_hypercube(dim),
+ quadrature,
+ face_no,
+ subface_no,
+ ref_case);
+}
+
+
+
+template <int dim>
+Quadrature<dim>
+QProjector<dim>::project_to_subface(
+ const ReferenceCell::Type reference_cell_type,
+ const SubQuadrature & quadrature,
+ const unsigned int face_no,
+ const unsigned int subface_no,
+ const RefinementCase<dim - 1> &ref_case)
+{
+ Assert(reference_cell_type == ReferenceCell::get_hypercube(dim),
+ ExcNotImplemented());
+ (void)reference_cell_type;
+
+ 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());
+}
+
+
+// explicit instantiations; note: we need them all for all dimensions
+template class QProjector<1>;
+template class QProjector<2>;
+template class QProjector<3>;
+
+DEAL_II_NAMESPACE_CLOSE
-//---------------------------------------------------------------------------
-
-
-
-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)
-{
- AssertIndexRange(child_no, 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
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
{
return get_data(flags,
mapping,
- QProjector<dim>::project_to_all_faces(quadrature),
+ QProjector<dim>::project_to_all_faces(
+ this->reference_cell_type(), quadrature),
output_data);
}
{
return get_data(flags,
mapping,
- QProjector<dim>::project_to_all_subfaces(quadrature),
+ QProjector<dim>::project_to_all_subfaces(
+ this->reference_cell_type(), quadrature),
output_data);
}
}
Quadrature<dim> faces =
- QProjector<dim>::project_to_all_faces(face_points);
+ QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
+ face_points);
for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
++current)
{
// child cell are evaluated
// in the quadrature points
// of a full face.
- Quadrature<dim> q_face = QProjector<dim>::project_to_face(q_base, face);
+ Quadrature<dim> q_face =
+ QProjector<dim>::project_to_face(this->reference_cell_type(),
+ q_base,
+ face);
// Store shape values, since the
// evaluation suffers if not
// ordered by point
// the coarse face are
// evaluated on the subface
// only.
- Quadrature<dim> q_sub =
- QProjector<dim>::project_to_subface(q_base, face, sub);
+ Quadrature<dim> q_sub = QProjector<dim>::project_to_subface(
+ this->reference_cell_type(), q_base, face, sub);
const unsigned int child = GeometryInfo<dim>::child_cell_on_face(
RefinementCase<dim>::isotropic_refinement, face, sub);
for (unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
++child)
{
- Quadrature<dim> q_sub = QProjector<dim>::project_to_child(q_cell, child);
+ Quadrature<dim> q_sub =
+ QProjector<dim>::project_to_child(this->reference_cell_type(),
+ q_cell,
+ child);
for (unsigned int k = 0; k < q_sub.size(); ++k)
for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
// TODO: Check what the face_orientation, face_flip and face_rotation
// have to be in 3D
unsigned int k = QProjector<dim>::DataSetDescriptor::face(
- face, false, false, false, n_face_points);
+ this->reference_cell_type(),
+ face,
+ false,
+ false,
+ false,
+ n_face_points);
for (unsigned int i = 0; i < boundary_weights_abf.size(1); ++i)
nodal_values[start_abf_dofs + i] +=
n_orient * boundary_weights_abf(k + fp, i) *
this->generalized_support_points.resize(npoints);
- Quadrature<dim> faces = QProjector<dim>::project_to_all_faces(face_points);
+ Quadrature<dim> faces =
+ QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
+ face_points);
for (unsigned int k = 0;
k < face_points.size() * GeometryInfo<dim>::faces_per_cell;
++k)
- this->generalized_support_points[k] =
- faces.point(k + QProjector<dim>::DataSetDescriptor::face(
- 0, true, false, false, this->n_dofs_per_face()));
+ this->generalized_support_points[k] = faces.point(
+ k + QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
+ 0,
+ true,
+ false,
+ false,
+ this->n_dofs_per_face()));
// Currently, for backward compatibility, we do not use moments, but
// point values on faces in 2D. In 3D, this is impossible, since the
// be done for the face orientation flag in 3D.
const Quadrature<dim> subface_quadrature =
subface == numbers::invalid_unsigned_int ?
- QProjector<dim>::project_to_face(quad_face_support, 0) :
- QProjector<dim>::project_to_subface(quad_face_support, 0, subface);
+ QProjector<dim>::project_to_face(this->reference_cell_type(),
+ quad_face_support,
+ 0) :
+ QProjector<dim>::project_to_subface(this->reference_cell_type(),
+ quad_face_support,
+ 0,
+ subface);
for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i)
{
const unsigned int n_boundary_points =
GeometryInfo<dim>::lines_per_cell * n_edge_points;
const Quadrature<dim> edge_quadrature =
- QProjector<dim>::project_to_all_faces(reference_edge_quadrature);
+ QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
+ reference_edge_quadrature);
this->generalized_face_support_points.resize(n_edge_points);
++line)
this->generalized_support_points[line * n_edge_points + q_point] =
edge_quadrature.point(QProjector<dim>::DataSetDescriptor::face(
- line, true, false, false, n_edge_points) +
+ this->reference_cell_type(),
+ line,
+ true,
+ false,
+ false,
+ n_edge_points) +
q_point);
for (unsigned int i = 0; i < order; ++i)
for (unsigned int q_point = 0; q_point < n_edge_points; ++q_point)
this->generalized_support_points[line * n_edge_points + q_point] =
edge_quadrature.point(QProjector<dim>::DataSetDescriptor::face(
- line, true, false, false, n_edge_points) +
+ this->reference_cell_type(),
+ line,
+ true,
+ false,
+ false,
+ n_edge_points) +
q_point);
}
}
const QGauss<1> reference_edge_quadrature(order + 1);
const unsigned int n_edge_points = reference_edge_quadrature.size();
const Quadrature<dim - 1> &edge_quadrature =
- QProjector<dim - 1>::project_to_all_faces(reference_edge_quadrature);
+ QProjector<dim - 1>::project_to_all_faces(ReferenceCell::get_hypercube(dim -
+ 1),
+ reference_edge_quadrature);
if (order > 0)
{
q_point] =
edge_quadrature.point(
QProjector<dim - 1>::DataSetDescriptor::face(
- line, true, false, false, n_edge_points) +
+ ReferenceCell::get_hypercube(dim - 1),
+ line,
+ true,
+ false,
+ false,
+ n_edge_points) +
q_point);
for (unsigned int i = 0; i < 2; ++i)
}
const Quadrature<dim> &face_quadrature =
- QProjector<dim>::project_to_all_faces(reference_face_quadrature);
+ QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
+ reference_face_quadrature);
for (const unsigned int face : GeometryInfo<dim>::face_indices())
for (unsigned int q_point = 0; q_point < n_face_points; ++q_point)
GeometryInfo<dim>::lines_per_cell *
n_edge_points] =
face_quadrature.point(QProjector<dim>::DataSetDescriptor::face(
- face, true, false, false, n_face_points) +
+ this->reference_cell_type(),
+ face,
+ true,
+ false,
+ false,
+ n_face_points) +
q_point);
}
q_point] =
edge_quadrature.point(
QProjector<dim - 1>::DataSetDescriptor::face(
- line, true, false, false, n_edge_points) +
+ ReferenceCell::get_hypercube(dim - 1),
+ line,
+ true,
+ false,
+ false,
+ n_edge_points) +
q_point);
for (unsigned int i = 0; i < 2; ++i)
// This will fill in the missing items in the InternalData
// (fe_internal/fe_data) which was not filled in by get_data.
fill_edge_values(cell,
- QProjector<dim>::project_to_all_faces(quadrature),
+ QProjector<dim>::project_to_all_faces(
+ this->reference_cell_type(), quadrature),
fe_data);
if (dim == 3 && this->degree > 1)
{
fill_face_values(cell,
- QProjector<dim>::project_to_all_faces(quadrature),
+ QProjector<dim>::project_to_all_faces(
+ this->reference_cell_type(), quadrature),
fe_data);
}
const UpdateFlags flags(fe_data.update_each);
const unsigned int n_q_points = quadrature.size();
- const typename QProjector<dim>::DataSetDescriptor offset =
- QProjector<dim>::DataSetDescriptor::face(face_no,
+ const auto offset =
+ QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
+ face_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
cell->face_rotation(face_no),
// compute on the face
const Quadrature<2> quadrature_on_face =
- QProjector<2>::project_to_face(quadrature, face_no);
+ QProjector<2>::project_to_face(this->reference_cell_type(),
+ quadrature,
+ face_no);
if (flags & update_values)
for (unsigned int i = 0; i < quadrature_on_face.size(); ++i)
get_linear_shape_coefficients(cell);
// compute on the subface
- const Quadrature<2> quadrature_on_subface =
- QProjector<2>::project_to_subface(quadrature, face_no, sub_no);
+ const Quadrature<2> quadrature_on_subface = QProjector<2>::project_to_subface(
+ this->reference_cell_type(), quadrature, face_no, sub_no);
if (flags & update_values)
for (unsigned int i = 0; i < quadrature_on_subface.size(); ++i)
// to take (all data sets for all
// faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset =
- QProjector<dim>::DataSetDescriptor::face(face_no,
+ const auto offset =
+ QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
+ face_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
cell->face_rotation(face_no),
// offset determines which data set
// to take (all data sets for all
// sub-faces are stored contiguously)
- const typename QProjector<dim>::DataSetDescriptor offset =
- QProjector<dim>::DataSetDescriptor::subface(face_no,
+ const auto offset =
+ QProjector<dim>::DataSetDescriptor::subface(this->reference_cell_type(),
+ face_no,
sub_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
// Add nodes of lines interior in the "mother-face"
// line 5: use line 9
- QProjector<dim - 1>::project_to_subface(qline, 0, 0, p_line);
+ QProjector<dim - 1>::project_to_subface(
+ ReferenceCell::get_hypercube(dim - 1), qline, 0, 0, p_line);
for (unsigned int i = 0; i < n; ++i)
constraint_points.push_back(p_line[i] + Point<dim - 1>(0.5, 0));
// line 6: use line 10
- QProjector<dim - 1>::project_to_subface(qline, 0, 1, p_line);
+ QProjector<dim - 1>::project_to_subface(
+ ReferenceCell::get_hypercube(dim - 1), qline, 0, 1, p_line);
for (unsigned int i = 0; i < n; ++i)
constraint_points.push_back(p_line[i] + Point<dim - 1>(0.5, 0));
// line 7: use line 13
- QProjector<dim - 1>::project_to_subface(qline, 2, 0, p_line);
+ QProjector<dim - 1>::project_to_subface(
+ ReferenceCell::get_hypercube(dim - 1), qline, 2, 0, p_line);
for (unsigned int i = 0; i < n; ++i)
constraint_points.push_back(p_line[i] + Point<dim - 1>(0, 0.5));
// line 8: use line 14
- QProjector<dim - 1>::project_to_subface(qline, 2, 1, p_line);
+ QProjector<dim - 1>::project_to_subface(
+ ReferenceCell::get_hypercube(dim - 1), qline, 2, 1, p_line);
for (unsigned int i = 0; i < n; ++i)
constraint_points.push_back(p_line[i] + Point<dim - 1>(0, 0.5));
subface < GeometryInfo<dim - 1>::max_children_per_face;
++subface)
{
- QProjector<dim - 1>::project_to_subface(qline,
- face,
- subface,
- p_line);
+ QProjector<dim - 1>::project_to_subface(
+ ReferenceCell::get_hypercube(dim - 1),
+ qline,
+ face,
+ subface,
+ p_line);
constraint_points.insert(constraint_points.end(),
p_line.begin(),
p_line.end());
// be done for the face orientation flag in 3D.
const Quadrature<dim> subface_quadrature =
subface == numbers::invalid_unsigned_int ?
- QProjector<dim>::project_to_face(quad_face_support, 0) :
- QProjector<dim>::project_to_subface(quad_face_support, 0, subface);
+ QProjector<dim>::project_to_face(this->reference_cell_type(),
+ quad_face_support,
+ 0) :
+ QProjector<dim>::project_to_subface(this->reference_cell_type(),
+ quad_face_support,
+ 0,
+ subface);
for (unsigned int i = 0; i < source_fe->n_dofs_per_face(); ++i)
{
const Point<dim> &p = subface_quadrature.point(i);
// ---------------------------------------------------------------------
-#include <deal.II/base/qprojector.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/template_constraints.h>
// ---------------------------------------------------------------------
-#include <deal.II/base/qprojector.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/template_constraints.h>
}
Quadrature<dim> faces =
- QProjector<dim>::project_to_all_faces(face_points);
+ QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
+ face_points);
for (; current < GeometryInfo<dim>::faces_per_cell * n_face_points;
++current)
{
// into the vector
this->generalized_support_points[current] =
faces.point(current + QProjector<dim>::DataSetDescriptor::face(
- 0, true, false, false, n_face_points));
+ this->reference_cell_type(),
+ 0,
+ true,
+ false,
+ false,
+ n_face_points));
}
}
// child cell are evaluated
// in the quadrature points
// of a full face.
- Quadrature<dim> q_face = QProjector<dim>::project_to_face(q_base, face);
+ Quadrature<dim> q_face =
+ QProjector<dim>::project_to_face(this->reference_cell_type(),
+ q_base,
+ face);
// Store shape values, since the
// evaluation suffers if not
// ordered by point
// the coarse face are
// evaluated on the subface
// only.
- Quadrature<dim> q_sub =
- QProjector<dim>::project_to_subface(q_base, face, sub);
+ Quadrature<dim> q_sub = QProjector<dim>::project_to_subface(
+ this->reference_cell_type(), q_base, face, sub);
const unsigned int child = GeometryInfo<dim>::child_cell_on_face(
RefinementCase<dim>::isotropic_refinement, face, sub);
for (unsigned int child = 0; child < GeometryInfo<dim>::max_children_per_cell;
++child)
{
- Quadrature<dim> q_sub = QProjector<dim>::project_to_child(q_cell, child);
+ Quadrature<dim> q_sub =
+ QProjector<dim>::project_to_child(this->reference_cell_type(),
+ q_cell,
+ child);
for (unsigned int k = 0; k < q_sub.size(); ++k)
for (unsigned int i_child = 0; i_child < this->n_dofs_per_cell();
for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
this->generalized_face_support_points[k] = face_points.point(k);
Quadrature<dim> faces =
- QProjector<dim>::project_to_all_faces(face_points);
+ QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
+ face_points);
for (unsigned int k = 0;
k < this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
++k)
- this->generalized_support_points[k] =
- faces.point(k + QProjector<dim>::DataSetDescriptor::face(
- 0, true, false, false, this->n_dofs_per_face()));
+ this->generalized_support_points[k] = faces.point(
+ k +
+ QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
+ 0,
+ true,
+ false,
+ false,
+ this->n_dofs_per_face()));
current = this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
}
// matrix by simply taking the
// value at the support points.
const Quadrature<dim> face_projection =
- QProjector<dim>::project_to_face(quad_face_support, 0);
+ QProjector<dim>::project_to_face(this->reference_cell_type(),
+ quad_face_support,
+ 0);
for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i)
{
// value at the support points.
const Quadrature<dim> subface_projection =
- QProjector<dim>::project_to_subface(quad_face_support, 0, subface);
+ QProjector<dim>::project_to_subface(this->reference_cell_type(),
+ quad_face_support,
+ 0,
+ subface);
for (unsigned int i = 0; i < source_fe.n_dofs_per_face(); ++i)
{
for (unsigned int k = 0; k < this->n_dofs_per_face(); ++k)
this->generalized_face_support_points[k] = face_points.point(k);
Quadrature<dim> faces =
- QProjector<dim>::project_to_all_faces(face_points);
+ QProjector<dim>::project_to_all_faces(this->reference_cell_type(),
+ face_points);
for (unsigned int k = 0;
k < this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
++k)
- this->generalized_support_points[k] =
- faces.point(k + QProjector<dim>::DataSetDescriptor::face(
- 0, true, false, false, this->n_dofs_per_face()));
+ this->generalized_support_points[k] = faces.point(
+ k +
+ QProjector<dim>::DataSetDescriptor::face(this->reference_cell_type(),
+ 0,
+ true,
+ false,
+ false,
+ this->n_dofs_per_face()));
current = this->n_dofs_per_face() * GeometryInfo<dim>::faces_per_cell;
}
// ---------------------------------------------------------------------
#include <deal.II/base/memory_consumption.h>
-#include <deal.II/base/qprojector.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/base/config.h>
-#include <deal.II/base/qprojector.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor_product_polynomials.h>
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(
- QProjector<dim>::project_to_all_faces(quadrature));
+ QProjector<dim>::project_to_all_faces(ReferenceCell::get_hypercube(dim),
+ quadrature));
auto &data = dynamic_cast<InternalData &>(*data_ptr);
// verify that we have computed the transitive hull of the required
const Quadrature<dim - 1> &quadrature) const
{
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
- std::make_unique<InternalData>(
- QProjector<dim>::project_to_all_subfaces(quadrature));
+ std::make_unique<InternalData>(QProjector<dim>::project_to_all_subfaces(
+ ReferenceCell::get_hypercube(dim), quadrature));
auto &data = dynamic_cast<InternalData &>(*data_ptr);
// verify that we have computed the transitive hull of the required
{
if (data.update_each & update_quadrature_points)
{
- const typename QProjector<dim>::DataSetDescriptor offset =
- QProjector<dim>::DataSetDescriptor::cell();
+ const auto offset = QProjector<dim>::DataSetDescriptor::cell();
transform_quadrature_points(cell, data, offset, quadrature_points);
}
if (data.update_each & update_quadrature_points)
{
- const typename QProjector<dim>::DataSetDescriptor offset =
- QProjector<dim>::DataSetDescriptor::face(face_no,
- cell->face_orientation(
- face_no),
- cell->face_flip(face_no),
- cell->face_rotation(face_no),
- quadrature_points.size());
+ const auto offset = QProjector<dim>::DataSetDescriptor::face(
+ ReferenceCell::get_hypercube(dim),
+ face_no,
+ cell->face_orientation(face_no),
+ cell->face_flip(face_no),
+ cell->face_rotation(face_no),
+ quadrature_points.size());
transform_quadrature_points(cell, data, offset, quadrature_points);
if (data.update_each & update_quadrature_points)
{
- const typename QProjector<dim>::DataSetDescriptor offset =
- QProjector<dim>::DataSetDescriptor::subface(
- face_no,
- sub_no,
- cell->face_orientation(face_no),
- cell->face_flip(face_no),
- cell->face_rotation(face_no),
- quadrature_points.size(),
- cell->subface_case(face_no));
+ const auto offset = QProjector<dim>::DataSetDescriptor::subface(
+ ReferenceCell::get_hypercube(dim),
+ face_no,
+ sub_no,
+ cell->face_orientation(face_no),
+ cell->face_flip(face_no),
+ cell->face_rotation(face_no),
+ quadrature_points.size(),
+ cell->subface_case(face_no));
transform_quadrature_points(cell, data, offset, quadrature_points);
}
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
auto & data = dynamic_cast<InternalData &>(*data_ptr);
- const Quadrature<dim> q(QProjector<dim>::project_to_all_faces(quadrature));
+ const Quadrature<dim> q(
+ QProjector<dim>::project_to_all_faces(ReferenceCell::get_hypercube(dim),
+ quadrature));
this->compute_face_data(update_flags, q, quadrature.size(), data);
return data_ptr;
std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
auto & data = dynamic_cast<InternalData &>(*data_ptr);
- const Quadrature<dim> q(QProjector<dim>::project_to_all_subfaces(quadrature));
+ const Quadrature<dim> q(
+ QProjector<dim>::project_to_all_subfaces(ReferenceCell::get_hypercube(dim),
+ quadrature));
this->compute_face_data(update_flags, q, quadrature.size(), data);
return data_ptr;
cell,
face_no,
numbers::invalid_unsigned_int,
- QProjector<dim>::DataSetDescriptor::face(face_no,
+ QProjector<dim>::DataSetDescriptor::face(ReferenceCell::get_hypercube(
+ dim),
+ face_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
cell->face_rotation(face_no),
cell,
face_no,
numbers::invalid_unsigned_int,
- QProjector<dim>::DataSetDescriptor::subface(face_no,
- subface_no,
- cell->face_orientation(
- face_no),
- cell->face_flip(face_no),
- cell->face_rotation(face_no),
- quadrature.size(),
- cell->subface_case(face_no)),
+ QProjector<dim>::DataSetDescriptor::subface(
+ ReferenceCell::get_hypercube(dim),
+ face_no,
+ subface_no,
+ cell->face_orientation(face_no),
+ cell->face_flip(face_no),
+ cell->face_rotation(face_no),
+ quadrature.size(),
+ cell->subface_case(face_no)),
quadrature,
data,
euler_dof_handler->get_fe(),
std::make_unique<InternalData>();
auto &data = dynamic_cast<InternalData &>(*data_ptr);
data.initialize_face(this->requires_update_flags(update_flags),
- QProjector<dim>::project_to_all_faces(quadrature),
+ QProjector<dim>::project_to_all_faces(
+ ReferenceCell::get_hypercube(dim), quadrature),
quadrature.size());
return data_ptr;
std::make_unique<InternalData>();
auto &data = dynamic_cast<InternalData &>(*data_ptr);
data.initialize_face(this->requires_update_flags(update_flags),
- QProjector<dim>::project_to_all_subfaces(quadrature),
+ QProjector<dim>::project_to_all_subfaces(
+ ReferenceCell::get_hypercube(dim), quadrature),
quadrature.size());
return data_ptr;
cell,
face_no,
numbers::invalid_unsigned_int,
- QProjector<dim>::DataSetDescriptor::face(face_no,
+ QProjector<dim>::DataSetDescriptor::face(ReferenceCell::get_hypercube(dim),
+ face_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
cell->face_rotation(face_no),
cell,
face_no,
subface_no,
- QProjector<dim>::DataSetDescriptor::subface(face_no,
+ QProjector<dim>::DataSetDescriptor::subface(ReferenceCell::get_hypercube(
+ dim),
+ face_no,
subface_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
#include <deal.II/base/derivative_form.h>
#include <deal.II/base/memory_consumption.h>
-#include <deal.II/base/qprojector.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/tensor_product_polynomials.h>
std::make_unique<InternalData>(polynomial_degree);
auto &data = dynamic_cast<InternalData &>(*data_ptr);
data.initialize_face(this->requires_update_flags(update_flags),
- QProjector<dim>::project_to_all_faces(quadrature),
+ QProjector<dim>::project_to_all_faces(
+ ReferenceCell::get_hypercube(dim), quadrature),
quadrature.size());
return data_ptr;
std::make_unique<InternalData>(polynomial_degree);
auto &data = dynamic_cast<InternalData &>(*data_ptr);
data.initialize_face(this->requires_update_flags(update_flags),
- QProjector<dim>::project_to_all_subfaces(quadrature),
+ QProjector<dim>::project_to_all_subfaces(
+ ReferenceCell::get_hypercube(dim), quadrature),
quadrature.size());
return data_ptr;
cell,
face_no,
numbers::invalid_unsigned_int,
- QProjector<dim>::DataSetDescriptor::face(face_no,
+ QProjector<dim>::DataSetDescriptor::face(ReferenceCell::get_hypercube(dim),
+ face_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
cell->face_rotation(face_no),
cell,
face_no,
subface_no,
- QProjector<dim>::DataSetDescriptor::subface(face_no,
+ QProjector<dim>::DataSetDescriptor::subface(ReferenceCell::get_hypercube(
+ dim),
+ face_no,
subface_no,
cell->face_orientation(face_no),
cell->face_flip(face_no),
std::vector<double> dummy_weights(n_points, 1. / n_points);
Quadrature<dim - 1> quadrature(boundary_points, dummy_weights);
- q_projector = QProjector<dim>::project_to_all_faces(quadrature);
+ q_projector =
+ QProjector<dim>::project_to_all_faces(ReferenceCell::Type::Quad,
+ quadrature);
}
for (const auto &cell : tria.active_cell_iterators())
projects / dealii.git / commitdiff