* Given an arbitrary quadrature formula, return one that chops the quadrature
* points above the hyper-plane defined by $\sum_i x_i = 1$. In other words,
* it extracts those quadrature points from the base formula that satisfy
- * $\sum_i (\mathbf x_q)_i \le 1$."
+ * $\sum_i (\mathbf x_q)_i \le 1+10^{-12}$."
*
* In general the resulting quadrature is not very useful, unless the
* quadrature you started from has been constructed specifically to integrate
- * over triangles or tetrahedra. This class ensures that the resulting
+ * over triangles or tetrahedra. This class only ensures that the resulting
* quadrature formula only has quadrature points in the reference simplex or on
* its boundary.
*
- * No transformation is applied to the weights. The weights referring to points
- * that live outside the reference simplex are simply discarded. If you chop a
- * quadrature formula that is good for both the lower and upper triangles
- * separately, then the weights should add up to 1/2 for each of the two
- * formulas.
+ * No transformation is applied to the weights, and the weights referring to points
+ * that live outside the reference simplex are simply discarded.
+ *
+ * The main use of this quadrature formula is not to chop tensor product
+ * quadratures. Ideally you should pass to this class a quadrature formula
+ * constructed directly using points and weights in the reference simplex,
+ * capable of integrating on triangles or tetrahedra.
+ *
+ * For finite elements based on quadrilaterals and hexahedra, a QSimplex
+ * quadrature formula is not very useful on its own. This class is typically
+ * used in conjuction with other classes, like QSplit, to patch the reference
+ * element using several QSimplex quadrature formulas.
+ *
+ * Such quadrature formulas are useful to integrate functions with
+ * singularities at certain points, or functions that present jumps along a
+ * co-dimension one surface inside the reference element, like in the extended
+ * finite element method.
*
* @author Luca Heltai, 2017.
*/
};
/**
- * A two dimensional simplex quadrature with zero Jacobian on vertex zero. The
- * quadrature is obtained through the following polar transformation:
+ * A quadrature that implements a polar transformation from a square to a
+ * triangle to integrate singularities in the origin of the reference simplex.
+ * The quadrature is obtained through the following polar transformation:
*
* \f[
* \begin{pmatrix}
*
* with determinant of the Jacobian equal to $J= \beta \hat \x^{2\beta-1}$.
* Such transformation maps the reference square \$[0,1]\times[0,1]$ to the
- * reference simplex, by collapsing the left \side of the square and squeezing
+ * reference simplex, by collapsing the left side of the square and squeezing
* quadrature points towards the orgin, and then shearing the resulting
* triangle to the reference one. This transformation shows good convergence
* properties when $\beta = 1$ with singularities of order $1/R$ in the origin,
*/
QDuffy(const Quadrature<1> &radial_quadrature,
const Quadrature<1> &angular_quadrature,
- const double &beta = 1.0);
+ const double beta = 1.0);
/**
- * Calls the above constructor with QGauss<1>(n) quadrature formulas for
+ * Call the above constructor with QGauss<1>(n) quadrature formulas for
* both the radial and angular quadratures.
*
* @param n
*/
- QDuffy(const unsigned int &n, const double &beta);
+ QDuffy(const unsigned int n,
+ const double beta);
};
/**
- * A quadrature to use when the cell should be split in subregions to integrate
- * using one or more base quadratures.
+ * A quadrature to use when the cell should be split into subregions to
+ * integrate using one or more base quadratures.
*
* @author Luca Heltai, 2017.
*/
* Construct a quadrature formula by splitting the reference hyper cube into
* the minimum number of simplices that have vertex zero coinciding with
* `split_point`, and patch together affine transformations of the `base`
- * quadrature.
+ * quadrature. The point `split_point` should be in the reference element,
+ * and an exception is thrown if this is not the case.
+ *
+ * In two dimensions, the resulting quadrature formula will be composed of
+ * two, three, or four triangular quadrature formulas if `split_point`
+ * coincides with one of the vertices, if it lies on one of the edges, or if
+ * it is internal to the reference element respectively.
+ *
+ * The same is true for the three dimensional case, with six, eight, ten, or
+ * twelve tetrahedral quadrature formulas if `split_point` coincides with one
+ * of the vertices, if it lies on one of the edges, on one of the faces, or
+ * if it is internal to the reference element respectively.
*
* The resulting quadrature can be used, for example, to integrate functions
* with integrable singularities at the split point, provided that you select
double r=0;
for (unsigned int d=0; d<dim; ++d)
r += quad.point(i)[d];
- if (r <= 1)
+ if (r <= 1+1e-10)
{
this->quadrature_points.push_back(quad.point(i));
this->weights.push_back(quad.weight(i));
B[d] = vertices[d+1]-vertices[0];
B = transpose(B);
- auto J = std::abs(determinant(B));
+ const double J = std::abs(determinant(B));
// if the determinant is zero, we return an empty quadrature
if (J < 1e-12)
const Quadrature<1> &angular_quadrature) :
QSimplex<2>(Quadrature<2>())
{
- QAnisotropic<2> base(radial_quadrature, angular_quadrature);
+ const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
this->quadrature_points.resize(base.size());
this->weights.resize(base.size());
for (unsigned int i=0; i<base.size(); ++i)
{
- const auto &q = base.point(i);
- const auto &w = base.weight(i);
+ const auto q = base.point(i);
+ const auto w = base.weight(i);
- const auto &xhat = q[0];
- const auto &yhat = q[1];
+ const auto xhat = q[0];
+ const auto yhat = q[1];
const double t = numbers::PI_2*yhat;
- const double &pi = numbers::PI;
+ const double pi = numbers::PI;
const double st = std::sin(t);
const double ct = std::cos(t);
const double r = xhat/(st+ct);
QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
const Quadrature<1> &angular_quadrature,
- const double &beta) :
+ const double beta) :
QSimplex<2>(Quadrature<2>())
{
- QAnisotropic<2> base(radial_quadrature, angular_quadrature);
+ const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
this->quadrature_points.resize(base.size());
this->weights.resize(base.size());
for (unsigned int i=0; i<base.size(); ++i)
{
- const auto &q = base.point(i);
- const auto &w = base.weight(i);
+ const auto q = base.point(i);
+ const auto w = base.weight(i);
- const auto &xhat = q[0];
- const auto &yhat = q[1];
+ const auto xhat = q[0];
+ const auto yhat = q[1];
const double x = std::pow(xhat, beta)*(1-yhat);
const double y = std::pow(xhat, beta)*yhat;
-QDuffy::QDuffy(const unsigned int &n, const double &beta)
+QDuffy::QDuffy(const unsigned int n, const double beta)
:QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
{}
QSplit<dim>::QSplit(const QSimplex<dim> &base,
const Point<dim> &split_point)
{
+ AssertThrow(GeometryInfo<dim>::is_inside_unit_cell(split_point, 1e-12),
+ ExcMessage("The split point should be inside the unit reference cell."));
+
std::array<Point<dim>, dim+1> vertices;
vertices[0] = split_point;
GeometryInfo<dim>::unit_cell_vertex(
GeometryInfo<dim>::face_to_cell_vertices(f, i+start)
);
- auto quad = base.compute_affine_transformation(vertices);
+ const auto quad = base.compute_affine_transformation(vertices);
if (quad.size())
{
this->quadrature_points.insert(
projects / dealii.git / commitdiff