<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
<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.
*/
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:
/**
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);
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))
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];
- // 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);
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);
*
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)
}
}
}
- // 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);
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)
{
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;
}
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)
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);
--- /dev/null
+//---------------------------- 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);
+}
+
+
+
--- /dev/null
+
+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))
--- /dev/null
+//------------------- 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;
+}
+
+
+
--- /dev/null
+
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::Integral = -1.0000
-//---------------------------- negative_jac_det.cc -------------------------
+//------------------- distorted_mapped_cells_01.cc ------------------------
// $Id$
// Version: $Name$
//
// 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();
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);
--- /dev/null
+
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::(typename Mapping<dim,spacedim>::ExcDistortedMappedCell(cell->center(), point))
+DEAL::Integral = 0.0757
+++ /dev/null
-
-DEAL::Integral = 1.0000, should be 1.