From: Luca Heltai Date: Thu, 21 Dec 2017 12:02:54 +0000 (+0100) Subject: Allow for negative volumes in QSimplex::compute_affine_transformation X-Git-Tag: v9.0.0-rc1~583^2~10 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a1c4a8b8bda2ecc739ad519de2ecc73608471e35;p=dealii.git Allow for negative volumes in QSimplex::compute_affine_transformation This is compensated by scaling the weights with the absolute values of the determinant of the Jacobian. --- diff --git a/include/deal.II/base/quadrature_lib.h b/include/deal.II/base/quadrature_lib.h index 03a1916961..5d213d901a 100644 --- a/include/deal.II/base/quadrature_lib.h +++ b/include/deal.II/base/quadrature_lib.h @@ -691,10 +691,11 @@ public: * \f] * where the matrix $B$ is given by $B_{ij} = v[j][i]-v[0][i]$. * - * The weights are scaled with the determinant of $B$. An assertion is thrown - * if the ordering of the vertices is such that the determinant of $B$ is - * negative (this means you attempted to construct a simplex with a negative - * area). + * The weights are scaled with the absolute value of the determinant of $B$, + * that is $J := |\text{det}(B)|$. If $J$ is zero, an empty quadrature is + * returned. This may happen, in two dimensions, if the three vertices are + * aligned, or in three dimensions if the four vertices are on the same + * plane. * * @param[in] vertices The vertices of the simplex you wish to integrate on * @return A quadrature object that can be used to integrate on the simplex diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc index 522f180988..9a3ae2c92a 100644 --- a/source/base/quadrature_lib.cc +++ b/source/base/quadrature_lib.cc @@ -1372,19 +1372,20 @@ template Quadrature QSimplex::compute_affine_transformation(const std::array, dim+1>& vertices) const { - std::vector > qp(this->size()); - std::vector w(this->size()); - unsigned int i=0; Tensor<2,dim> B; for (unsigned int d=0; d0, ExcMessage("The provided vertices generate an inverted Simplex." - " Make sure the vertices are oriented correctly.")) + // if the determinant is zero, we return an empty quadrature + if (J < 1e-12) + return Quadrature(); + + std::vector > qp(this->size()); + std::vector w(this->size()); for (unsigned int i=0; isize(); ++i) {