]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Check that the Jacobian determinant is always positive.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Jan 2013 08:08:47 +0000 (08:08 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 18 Jan 2013 08:08:47 +0000 (08:08 +0000)
git-svn-id: https://svn.dealii.org/trunk@28119 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/changes.h
deal.II/include/deal.II/fe/mapping.h
deal.II/source/fe/mapping_q1.cc
tests/bits/distorted_cells_07.cc [new file with mode: 0644]
tests/bits/distorted_cells_07/cmp/generic [new file with mode: 0644]
tests/bits/distorted_mapped_cells_01.cc [new file with mode: 0644]
tests/bits/distorted_mapped_cells_01/cmp/generic [new file with mode: 0644]
tests/bits/distorted_mapped_cells_02.cc [moved from tests/fe/negative_jac_det.cc with 52% similarity]
tests/bits/distorted_mapped_cells_02/cmp/generic [new file with mode: 0644]
tests/fe/negative_jac_det/cmp/generic [deleted file]

index 8406fdc8d38d3f746c022a4f6d8537d45088d0d4..90ed72b4b86a29250466db51e699ce874986b3a4 100644 (file)
@@ -177,6 +177,14 @@ DoFHandler, in particular removal of specializations.
 <h3>Specific improvements</h3>
 
 <ol>
+<li> Fixed: The FEValues machinery silently accepted the case when the
+mapped cell (or the cell geometry) were distorted. An assertion has been
+added to the computation of the Jacobian determinants for the volume
+element that aborts whenever the Jacobian determinant in a quadrature
+point becomes too small or negative.
+<br>
+(Wolfgang Bangerth, 2013/01/17)
+
 <li> Fixed: Various variants of the TrilinosWrappers::SparseMatrix::reinit
 functions take a parameter <code>drop_tolerance</code> that allows to remove
 small values from the matrix and replace them by zero instead. This was not
@@ -252,6 +260,13 @@ is now fixed.
 <br>
 (Wolfgang Bangerth, 2012/11/20)
 
+<li> Improved: The inner product, norm, and mean value computation
+of deal.II's own Vector class are now parallelized by threads.
+The parallelization does not change the order in which the additions
+take place, ensuring exact reproducibility from one run to the next.
+<br>
+(Martin Kronbichler, 2012/11/18)
+
 <li> New: The TrilinosWrappers::PreconditionBase class now has
 a function TrilinosWrappers::PreconditionBase::Tvmult that
 allows applying the transpose preconditioner.
index 3a5dc6077707aad3792e623761a82daf4ac940a3..bbae280d1d6e58e9785a5412b72a51af5315e977 100644 (file)
@@ -718,6 +718,20 @@ public:
    */
   DeclException0(ExcTransformationFailed);
 
+  /**
+   * deal.II assumes the Jacobian determinant to be positive. When the cell
+   * geometry is distorted under the image of the mapping, the mapping becomes
+   * invalid and this exception is thrown.
+   *
+   * @ingroup Exceptions
+   */
+  DeclException2 (ExcDistortedMappedCell,
+                  Point<spacedim>, int,
+                  << "The image of the mapping applied to cell with center ["
+                  << arg1 << "] is distorted. The cell geometry or the "
+                  << "mapping are invalid and yield non-positive volume "
+                  << "fractions in quadrature point " << arg2 << ".");
+
 private:
 
   /**
index 211ad2c8b91a7a887c97ea07fc46a42da045849c..336eedd7403a26b2be53d9582730a5d175697d83 100644 (file)
@@ -739,8 +739,7 @@ MappingQ1<dim,spacedim>::fill_fe_values (
   std::vector<Point<spacedim> >                             &normal_vectors,
   CellSimilarity::Similarity                                &cell_similarity) const
 {
-  // ensure that the following cast
-  // is really correct:
+  // ensure that the following static_cast is really correct:
   Assert (dynamic_cast<InternalData *>(&mapping_data) != 0,
           ExcInternalError());
   InternalData &data = static_cast<InternalData &>(mapping_data);
@@ -754,10 +753,8 @@ MappingQ1<dim,spacedim>::fill_fe_values (
   const UpdateFlags update_flags(data.current_update_flags());
   const std::vector<double> &weights=q.get_weights();
 
-  // Multiply quadrature weights by
-  // Jacobian determinants or the area
-  // element g=sqrt(DX^t DX) in case of
-  // codim > 0
+  // Multiply quadrature weights by absolute value of Jacobian determinants or
+  // the area element g=sqrt(DX^t DX) in case of codim > 0
 
   if (update_flags & (update_normal_vectors
                       | update_JxW_values))
@@ -775,18 +772,22 @@ MappingQ1<dim,spacedim>::fill_fe_values (
 
             if (dim == spacedim)
               {
-                JxW_values[point]
-                  = data.contravariant[point].determinant() * weights[point];
+                const double det = data.contravariant[point].determinant();
+
+                // check for distorted cells.
+
+                // TODO: this allows for anisotropies of up to 1e6 in 3D and
+                // 1e12 in 2D. might want to find a finer
+                // (dimension-independent) criterion
+                Assert (det > 1e-12*Utilities::fixed_power<dim>(cell->diameter()/
+                                                                std::sqrt(dim)),
+                        (typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point)));
+
+                JxW_values[point] = weights[point] * det;
               }
-            // if dim==spacedim,
-            // then there is no
-            // cell normal to
-            // compute. since this
-            // is for FEValues (and
-            // not FEFaceValues),
-            // there are also no
-            // face normals to
-            // compute
+            // if dim==spacedim, then there is no cell normal to
+            // compute. since this is for FEValues (and not FEFaceValues),
+            // there are also no face normals to compute
             else //codim>0 case
               {
                 Tensor<1, spacedim> DX_t [dim];
@@ -836,8 +837,7 @@ MappingQ1<dim,spacedim>::fill_fe_values (
 
 
 
-  // copy values from InternalData to
-  // vector given by reference
+  // copy values from InternalData to vector given by reference
   if (update_flags & update_jacobians)
     {
       AssertDimension (jacobians.size(), n_q_points);
@@ -846,10 +846,8 @@ MappingQ1<dim,spacedim>::fill_fe_values (
           jacobians[point] = data.contravariant[point];
     }
 
-  // calculate values of the
-  // derivatives of the Jacobians. do
-  // it here, since we only do it for
-  // cells, not faces.
+  // calculate values of the derivatives of the Jacobians. do it here, since
+  // we only do it for cells, not faces.
   if (update_flags & update_jacobian_grads)
     {
       AssertDimension (jacobian_grads.size(), n_q_points);
@@ -883,9 +881,8 @@ MappingQ1<dim,spacedim>::fill_fe_values (
                           *
                           data.mapping_support_points[k][i]);
 
-              // never touch any data for j=dim in case
-              // dim<spacedim, so it will always be zero as
-              // it was initialized
+              // never touch any data for j=dim in case dim<spacedim, so it
+              // will always be zero as it was initialized
               for (unsigned int i=0; i<spacedim; ++i)
                 for (unsigned int j=0; j<dim; ++j)
                   for (unsigned int l=0; l<dim; ++l)
@@ -893,8 +890,7 @@ MappingQ1<dim,spacedim>::fill_fe_values (
             }
         }
     }
-  // copy values from InternalData to vector
-  // given by reference
+  // copy values from InternalData to vector given by reference
   if (update_flags & update_inverse_jacobians)
     {
       AssertDimension (inverse_jacobians.size(), n_q_points);
@@ -937,10 +933,8 @@ namespace internal
           if (update_flags & update_JxW_values)
             AssertDimension (JxW_values.size(), n_q_points);
 
-          // map the unit tangentials to the
-          // real cell. checking for d!=dim-1
-          // eliminates compiler warnings
-          // regarding unsigned int expressions <
+          // map the unit tangentials to the real cell. checking for d!=dim-1
+          // eliminates compiler warnings regarding unsigned int expressions <
           // 0.
           for (unsigned int d=0; d!=dim-1; ++d)
             {
@@ -957,30 +951,18 @@ namespace internal
                                  mapping_contravariant);
             }
 
-          // if dim==spacedim, we can use the
-          // unit tangentials to compute the
-          // boundary form by simply taking
-          // the cross product
+          // if dim==spacedim, we can use the unit tangentials to compute the
+          // boundary form by simply taking the cross product
           if (dim == spacedim)
             {
               for (unsigned int i=0; i<n_q_points; ++i)
                 switch (dim)
                   {
                   case 1:
-                    // in 1d, we don't
-                    // have access to
-                    // any of the
-                    // data.aux fields
-                    // (because it has
-                    // only dim-1
-                    // components), but
-                    // we can still
-                    // compute the
-                    // boundary form
-                    // simply by simply
-                    // looking at the
-                    // number of the
-                    // face
+                    // in 1d, we don't have access to any of the data.aux
+                    // fields (because it has only dim-1 components), but we
+                    // can still compute the boundary form simply by simply
+                    // looking at the number of the face
                     boundary_forms[i][0] = (face_no == 0 ?
                                             -1 : +1);
                     break;
@@ -996,16 +978,12 @@ namespace internal
             }
           else //(dim < spacedim)
             {
-              // in the codim-one case, the
-              // boundary form results from
-              // the cross product of all the
-              // face tangential vectors and
-              // the cell normal vector
+              // in the codim-one case, the boundary form results from the
+              // cross product of all the face tangential vectors and the cell
+              // normal vector
               //
-              // to compute the cell normal,
-              // use the same method used in
-              // fill_fe_values for cells
-              // above
+              // to compute the cell normal, use the same method used in
+              // fill_fe_values for cells above
               AssertDimension (data.contravariant.size(), n_q_points);
 
               for (unsigned int point=0; point<n_q_points; ++point)
@@ -1027,10 +1005,8 @@ namespace internal
                       cross_product(cell_normal,DX_t[0],DX_t[1]);
                       cell_normal /= cell_normal.norm();
 
-                      // then compute the face
-                      // normal from the face
-                      // tangent and the cell
-                      // normal:
+                      // then compute the face normal from the face tangent
+                      // and the cell normal:
                       cross_product (boundary_forms[point],
                                      data.aux[0][point], cell_normal);
 
diff --git a/tests/bits/distorted_cells_07.cc b/tests/bits/distorted_cells_07.cc
new file mode 100644 (file)
index 0000000..f248665
--- /dev/null
@@ -0,0 +1,98 @@
+//----------------------------  distorted_cells_07.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2003, 2004, 2005, 2009, 2010, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  distorted_cells_07.cc  ---------------------------
+
+
+// check that the mapping throws an exception for the test case in distorted_cells_01.cc
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/tria_iterator.h>
+#include <deal.II/grid/grid_reordering.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
+
+#include <fstream>
+
+
+// create a (i) pinched cell (where two vertices coincide), or (ii)
+// twisted cell (where two vertices are swapped)
+template <int dim>
+void check (const unsigned int testcase)
+{
+  std::vector<Point<dim> > vertices;
+  for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
+    vertices.push_back (GeometryInfo<dim>::unit_cell_vertex(v));
+
+  switch (testcase)
+    {
+      case 2:
+           deallog << "Twisted cell in " << dim << "d" << std::endl;
+            std::swap (vertices[1], vertices[0]);
+           break;
+      default:
+           Assert (false, ExcNotImplemented());
+    }
+
+
+  std::vector<CellData<dim> > cells;
+  {
+    CellData<dim> cell;
+    for (unsigned int j=0; j<GeometryInfo<dim>::vertices_per_cell; ++j)
+      cell.vertices[j]   = j;
+    cells.push_back (cell);
+  }
+
+  Triangulation<dim> coarse_grid (Triangulation<dim>::none, true);
+
+  bool flag = false;
+  try
+    {
+      coarse_grid.create_triangulation (vertices, cells, SubCellData());
+    }
+  catch (typename Triangulation<dim>::DistortedCellList &dcv)
+    {
+      flag = true;
+    }
+
+  Assert (flag == true, ExcInternalError());
+
+  // now build an FEValues object and compute quadrature points on that cell
+  FE_Nothing<dim> dummy;
+  QGauss<dim> quadrature(2);
+  FEValues<dim> fe_values(dummy, quadrature, update_JxW_values);
+  // should throw an assertion
+  fe_values.reinit(coarse_grid.begin());
+}
+
+
+int main ()
+{
+  deal_II_exceptions::disable_abort_on_exception();
+  std::ofstream logfile("distorted_cells_07/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  // only twisted cells for FEValues (pinched cells are OK on Gauss points)
+  check<1> (2);
+  check<2> (2);
+  check<3> (2);
+}
+
+
+
diff --git a/tests/bits/distorted_cells_07/cmp/generic b/tests/bits/distorted_cells_07/cmp/generic
new file mode 100644 (file)
index 0000000..69c4e0c
--- /dev/null
@@ -0,0 +1,10 @@
+
+DEAL::Twisted cell in 1d
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::Twisted cell in 2d
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::Twisted cell in 3d
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
diff --git a/tests/bits/distorted_mapped_cells_01.cc b/tests/bits/distorted_mapped_cells_01.cc
new file mode 100644 (file)
index 0000000..dde1385
--- /dev/null
@@ -0,0 +1,87 @@
+//-------------------  distorted_mapped_cells_01.cc  ------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-------------------  distorted_mapped_cells_01.cc  ------------------------
+
+// Test that an assertion is thrown when MappingQEulerian produces a cell with
+// negative volume
+
+#include "../tests.h"
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+#include <deal.II/grid/tria.h>
+
+
+void test()
+{
+  const int dim = 2;
+  // create a dummy triangulation with no extension and set the geometry
+  // through MappingQEulerian
+  Triangulation<dim> tria;
+  std::vector<Point<dim> > points (4);
+
+  std::vector<CellData<dim> > cells (1);
+  cells[0].vertices[0] = 0;
+  cells[0].vertices[1] = 1;
+  cells[0].vertices[2] = 2;
+  cells[0].vertices[3] = 3;
+  cells[0].material_id = 0;
+
+  tria.create_triangulation(points, cells, SubCellData());
+
+  FE_Q<dim> fe(1);
+  FESystem<dim> fe_sys(fe, dim);
+  DoFHandler<dim> dof_h(tria);
+  dof_h.distribute_dofs(fe_sys);
+  Vector<double> displacements(dof_h.n_dofs());
+  displacements(2) = -1.;
+  displacements(5) = 1.;
+  displacements(6) = -1.;
+  displacements(7) = 1.;
+
+  // this gives a Cartesian cell but in non-standard orientation (x-coordinate
+  // is gone through backwards)
+  MappingQEulerian<dim> mapping(1,displacements, dof_h);
+  QGauss<dim> quad(1);
+  FEValues<dim> fe_val (mapping, fe, quad, update_JxW_values);
+  double integral = 0.;
+  /*typename*/ Triangulation<dim>::active_cell_iterator
+    cell = tria.begin_active(), endc = tria.end();
+  for ( ; cell != endc; ++cell)
+    {
+      fe_val.reinit (cell);
+      for (unsigned int q=0; q<quad.size(); ++q)
+        integral += fe_val.JxW(q);
+    }
+  deallog << "Integral = " << integral << std::endl;
+}
+
+
+int
+main()
+{
+  deal_II_exceptions::disable_abort_on_exception();
+  std::ofstream logfile ("distorted_mapped_cells_01/output");
+  deallog << std::setprecision(4) << std::fixed;
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  test();
+
+  return 0;
+}
+
+
+
diff --git a/tests/bits/distorted_mapped_cells_01/cmp/generic b/tests/bits/distorted_mapped_cells_01/cmp/generic
new file mode 100644 (file)
index 0000000..c8640ad
--- /dev/null
@@ -0,0 +1,3 @@
+
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::Integral = -1.0000
similarity index 52%
rename from tests/fe/negative_jac_det.cc
rename to tests/bits/distorted_mapped_cells_02.cc
index 766bf79c84d9308f48633aff5e14a71df6c6b7df..2cf6374f6f02ca9ee51e4d72e72ee57d558b6a31 100644 (file)
@@ -1,4 +1,4 @@
-//----------------------------  negative_jac_det.cc  -------------------------
+//-------------------  distorted_mapped_cells_01.cc  ------------------------
 //    $Id$
 //    Version: $Name$
 //
@@ -9,41 +9,35 @@
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  negative_jac_det.cc  -------------------------
+//-------------------  distorted_mapped_cells_01.cc  ------------------------
 
-// Test that the computation of the Jacobian determinant is correct also when
-// the mesh is in a non-standard orientation
+// Test that an assertion is thrown when MappingQ distorts the geometry too
+// much on a thin hyper shell
 
 #include "../tests.h"
 #include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/mapping_q.h>
 #include <deal.II/fe/fe_nothing.h>
 #include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/mapping_q1.h>
 #include <deal.II/grid/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria_boundary_lib.h>
 
 
 void test()
 {
   const int dim = 2;
-  Triangulation<dim> tria;
-  std::vector<Point<dim> > points (4);
-  points[0] = (Point<dim> (0, 0));
-  points[1] = (Point<dim> (0, 1));
-  points[2] = (Point<dim> (1, 0));
-  points[3] = (Point<dim> (1, 1));
+  Triangulation<dim> tria(Triangulation<dim>::none,true);
 
-  std::vector<CellData<dim> > cells (1);
-  cells[0].vertices[0] = 0;
-  cells[0].vertices[1] = 1;
-  cells[0].vertices[2] = 2;
-  cells[0].vertices[3] = 3;
-  cells[0].material_id = 0;
-
-  tria.create_triangulation(points, cells, SubCellData());
+  // shell is so narrow that MappingQ(2) distorts the cell
+  GridGenerator::quarter_hyper_shell (tria, Point<dim>(), 0.95, 1, 1);
+  static HyperShellBoundary<dim> boundary;
+  tria.set_boundary (0, boundary);
 
   FE_Nothing<dim> dummy;
+  MappingQ<dim> mapping(2);
   QGauss<dim> quad(2);
-  FEValues<dim> fe_val (dummy, quad, update_JxW_values);
+  FEValues<dim> fe_val (mapping, dummy, quad, update_JxW_values);
   double integral = 0.;
   /*typename*/ Triangulation<dim>::active_cell_iterator
     cell = tria.begin_active(), endc = tria.end();
@@ -53,14 +47,15 @@ void test()
       for (unsigned int q=0; q<quad.size(); ++q)
         integral += fe_val.JxW(q);
     }
-  deallog << "Integral = " << integral << ", should be 1." << std::endl;
+  deallog << "Integral = " << integral << std::endl;
 }
 
 
 int
 main()
 {
-  std::ofstream logfile ("negative_jac_det/output");
+  deal_II_exceptions::disable_abort_on_exception();
+  std::ofstream logfile ("distorted_mapped_cells_02/output");
   deallog << std::setprecision(4) << std::fixed;
   deallog.attach(logfile);
   deallog.depth_console(0);
diff --git a/tests/bits/distorted_mapped_cells_02/cmp/generic b/tests/bits/distorted_mapped_cells_02/cmp/generic
new file mode 100644 (file)
index 0000000..3c0a207
--- /dev/null
@@ -0,0 +1,4 @@
+
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::Integral = 0.0757
diff --git a/tests/fe/negative_jac_det/cmp/generic b/tests/fe/negative_jac_det/cmp/generic
deleted file mode 100644 (file)
index 264b0b7..0000000
+++ /dev/null
@@ -1,2 +0,0 @@
-
-DEAL::Integral = 1.0000, should be 1.

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.