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 03a1916961e2e14fb2c937fdfa84d29910eac85b..5d213d901a06a4bab2fdf8c6bfeb3f45fde26a26 100644 (file)
--- 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 522f180988fc84b7fc0d21bf3def4edab0df505c..9a3ae2c92a7776349c5f8316a23b1cb75e94b31b 100644 (file)
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -1372,19 +1372,20 @@ template<int dim>
 Quadrature<dim>
 QSimplex<dim>::compute_affine_transformation(const std::array<Point<dim>, dim+1>& vertices) const
 {
-  std::vector<Point<dim> > qp(this->size());
-  std::vector<double> w(this->size());
-
   unsigned int i=0;
   Tensor<2,dim> B;
   for (unsigned int d=0; d<dim; ++d)
     B[d] = vertices[d+1]-vertices[0];
 
   B = transpose(B);
-  auto J = determinant(B);
+  auto J = std::abs(determinant(B));
 
-  AssertThrow(J>0, 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<dim>();
+
+  std::vector<Point<dim> > qp(this->size());
+  std::vector<double> w(this->size());
 
   for (unsigned int i=0; i<this->size(); ++i)
     {