From d7b20d1f194a2e1051a6f5fe1a354c2b476ad109 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 18 Jan 2013 08:08:47 +0000 Subject: [PATCH] Check that the Jacobian determinant is always positive. git-svn-id: https://svn.dealii.org/trunk@28119 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 15 +++ deal.II/include/deal.II/fe/mapping.h | 14 +++ deal.II/source/fe/mapping_q1.cc | 102 +++++++----------- tests/bits/distorted_cells_07.cc | 98 +++++++++++++++++ tests/bits/distorted_cells_07/cmp/generic | 10 ++ tests/bits/distorted_mapped_cells_01.cc | 87 +++++++++++++++ .../distorted_mapped_cells_01/cmp/generic | 3 + .../distorted_mapped_cells_02.cc} | 39 +++---- .../distorted_mapped_cells_02/cmp/generic | 4 + tests/fe/negative_jac_det/cmp/generic | 2 - 10 files changed, 287 insertions(+), 87 deletions(-) create mode 100644 tests/bits/distorted_cells_07.cc create mode 100644 tests/bits/distorted_cells_07/cmp/generic create mode 100644 tests/bits/distorted_mapped_cells_01.cc create mode 100644 tests/bits/distorted_mapped_cells_01/cmp/generic rename tests/{fe/negative_jac_det.cc => bits/distorted_mapped_cells_02.cc} (52%) create mode 100644 tests/bits/distorted_mapped_cells_02/cmp/generic delete mode 100644 tests/fe/negative_jac_det/cmp/generic diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 8406fdc8d3..90ed72b4b8 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -177,6 +177,14 @@ DoFHandler, in particular removal of specializations.

Specific improvements

    +
  1. 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. +
    +(Wolfgang Bangerth, 2013/01/17) +
  2. Fixed: Various variants of the TrilinosWrappers::SparseMatrix::reinit functions take a parameter drop_tolerance 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.
    (Wolfgang Bangerth, 2012/11/20) +
  3. 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. +
    +(Martin Kronbichler, 2012/11/18) +
  4. New: The TrilinosWrappers::PreconditionBase class now has a function TrilinosWrappers::PreconditionBase::Tvmult that allows applying the transpose preconditioner. diff --git a/deal.II/include/deal.II/fe/mapping.h b/deal.II/include/deal.II/fe/mapping.h index 3a5dc60777..bbae280d1d 100644 --- a/deal.II/include/deal.II/fe/mapping.h +++ b/deal.II/include/deal.II/fe/mapping.h @@ -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, 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: /** diff --git a/deal.II/source/fe/mapping_q1.cc b/deal.II/source/fe/mapping_q1.cc index 211ad2c8b9..336eedd740 100644 --- a/deal.II/source/fe/mapping_q1.cc +++ b/deal.II/source/fe/mapping_q1.cc @@ -739,8 +739,7 @@ MappingQ1::fill_fe_values ( std::vector > &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(&mapping_data) != 0, ExcInternalError()); InternalData &data = static_cast(mapping_data); @@ -754,10 +753,8 @@ MappingQ1::fill_fe_values ( const UpdateFlags update_flags(data.current_update_flags()); const std::vector &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::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(cell->diameter()/ + std::sqrt(dim)), + (typename Mapping::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::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::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::fill_fe_values ( * data.mapping_support_points[k][i]); - // never touch any data for j=dim in case - // dim::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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +// create a (i) pinched cell (where two vertices coincide), or (ii) +// twisted cell (where two vertices are swapped) +template +void check (const unsigned int testcase) +{ + std::vector > vertices; + for (unsigned int v=0; v::vertices_per_cell; ++v) + vertices.push_back (GeometryInfo::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 > cells; + { + CellData cell; + for (unsigned int j=0; j::vertices_per_cell; ++j) + cell.vertices[j] = j; + cells.push_back (cell); + } + + Triangulation coarse_grid (Triangulation::none, true); + + bool flag = false; + try + { + coarse_grid.create_triangulation (vertices, cells, SubCellData()); + } + catch (typename Triangulation::DistortedCellList &dcv) + { + flag = true; + } + + Assert (flag == true, ExcInternalError()); + + // now build an FEValues object and compute quadrature points on that cell + FE_Nothing dummy; + QGauss quadrature(2); + FEValues 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 index 0000000000..69c4e0c314 --- /dev/null +++ b/tests/bits/distorted_cells_07/cmp/generic @@ -0,0 +1,10 @@ + +DEAL::Twisted cell in 1d +DEAL::(typename Mapping::ExcDistortedMappedCell(cell->center(), point)) +DEAL::(typename Mapping::ExcDistortedMappedCell(cell->center(), point)) +DEAL::Twisted cell in 2d +DEAL::(typename Mapping::ExcDistortedMappedCell(cell->center(), point)) +DEAL::(typename Mapping::ExcDistortedMappedCell(cell->center(), point)) +DEAL::Twisted cell in 3d +DEAL::(typename Mapping::ExcDistortedMappedCell(cell->center(), point)) +DEAL::(typename Mapping::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 index 0000000000..dde1385ddc --- /dev/null +++ b/tests/bits/distorted_mapped_cells_01.cc @@ -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 +#include +#include +#include +#include +#include + + +void test() +{ + const int dim = 2; + // create a dummy triangulation with no extension and set the geometry + // through MappingQEulerian + Triangulation tria; + std::vector > points (4); + + std::vector > 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 fe(1); + FESystem fe_sys(fe, dim); + DoFHandler dof_h(tria); + dof_h.distribute_dofs(fe_sys); + Vector 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 mapping(1,displacements, dof_h); + QGauss quad(1); + FEValues fe_val (mapping, fe, quad, update_JxW_values); + double integral = 0.; + /*typename*/ Triangulation::active_cell_iterator + cell = tria.begin_active(), endc = tria.end(); + for ( ; cell != endc; ++cell) + { + fe_val.reinit (cell); + for (unsigned int q=0; q::ExcDistortedMappedCell(cell->center(), point)) +DEAL::Integral = -1.0000 diff --git a/tests/fe/negative_jac_det.cc b/tests/bits/distorted_mapped_cells_02.cc similarity index 52% rename from tests/fe/negative_jac_det.cc rename to tests/bits/distorted_mapped_cells_02.cc index 766bf79c84..2cf6374f6f 100644 --- a/tests/fe/negative_jac_det.cc +++ b/tests/bits/distorted_mapped_cells_02.cc @@ -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 +#include #include #include -#include #include +#include +#include void test() { const int dim = 2; - Triangulation tria; - std::vector > points (4); - points[0] = (Point (0, 0)); - points[1] = (Point (0, 1)); - points[2] = (Point (1, 0)); - points[3] = (Point (1, 1)); + Triangulation tria(Triangulation::none,true); - std::vector > 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(), 0.95, 1, 1); + static HyperShellBoundary boundary; + tria.set_boundary (0, boundary); FE_Nothing dummy; + MappingQ mapping(2); QGauss quad(2); - FEValues fe_val (dummy, quad, update_JxW_values); + FEValues fe_val (mapping, dummy, quad, update_JxW_values); double integral = 0.; /*typename*/ Triangulation::active_cell_iterator cell = tria.begin_active(), endc = tria.end(); @@ -53,14 +47,15 @@ void test() for (unsigned int q=0; q::ExcDistortedMappedCell(cell->center(), point)) +DEAL::(typename Mapping::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 index 264b0b73cc..0000000000 --- a/tests/fe/negative_jac_det/cmp/generic +++ /dev/null @@ -1,2 +0,0 @@ - -DEAL::Integral = 1.0000, should be 1. -- 2.39.5