From: Simon Sticko Date: Tue, 17 May 2022 06:41:09 +0000 (+0200) Subject: Assert that the incoming cell is Cartesian in MappingCartesian X-Git-Tag: v9.4.0-rc1~205^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F13748%2Fhead;p=dealii.git Assert that the incoming cell is Cartesian in MappingCartesian --- diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc index 84cc8a3ddc..714553d7d3 100644 --- a/source/fe/mapping_cartesian.cc +++ b/source/fe/mapping_cartesian.cc @@ -37,6 +37,34 @@ DEAL_II_NAMESPACE_OPEN +DeclExceptionMsg( + ExcCellNotCartesian, + "You are using MappingCartesian, but the incoming cell is not Cartesian."); + + + +/** + * Return whether the incoming cell is of Cartesian shape. This is determined by + * checking if the smallest BoundingBox that encloses the cell has the same + * vertices as the cell itself. + */ +template +bool +is_cartesian(const CellType &cell) +{ + if (!cell->reference_cell().is_hyper_cube()) + return false; + + const double tolerance = 1e-14 * cell->diameter(); + const auto bounding_box = cell->bounding_box(); + + for (const unsigned int v : cell->vertex_indices()) + if (cell->vertex(v).distance(bounding_box.vertex(v)) > tolerance) + return false; + + return true; +} + template @@ -455,6 +483,8 @@ MappingCartesian::fill_fe_values( internal::FEValuesImplementation::MappingRelatedData &output_data) const { + Assert(is_cartesian(cell), ExcCellNotCartesian()); + // convert data object to internal data for this class. fails with // an exception if that is not possible Assert(dynamic_cast(&internal_data) != nullptr, @@ -504,6 +534,8 @@ MappingCartesian::fill_mapping_data_for_generic_points( if (update_flags == update_default) return; + Assert(is_cartesian(cell), ExcCellNotCartesian()); + Assert(update_flags & update_inverse_jacobians || update_flags & update_jacobians || update_flags & update_quadrature_points, @@ -538,6 +570,7 @@ MappingCartesian::fill_fe_face_values( internal::FEValuesImplementation::MappingRelatedData &output_data) const { + Assert(is_cartesian(cell), ExcCellNotCartesian()); AssertDimension(quadrature.size(), 1); // convert data object to internal @@ -591,6 +624,8 @@ MappingCartesian::fill_fe_subface_values( internal::FEValuesImplementation::MappingRelatedData &output_data) const { + Assert(is_cartesian(cell), ExcCellNotCartesian()); + // convert data object to internal data for this class. fails with // an exception if that is not possible Assert(dynamic_cast(&internal_data) != nullptr, @@ -647,6 +682,7 @@ MappingCartesian::fill_fe_immersed_surface_values( &output_data) const { AssertDimension(dim, spacedim); + Assert(is_cartesian(cell), ExcCellNotCartesian()); // Convert data object to internal data for this class. Fails with an // exception if that is not possible. @@ -1113,6 +1149,8 @@ MappingCartesian::transform_unit_to_real_cell( const typename Triangulation::cell_iterator &cell, const Point & p) const { + Assert(is_cartesian(cell), ExcCellNotCartesian()); + Tensor<1, dim> length; const Point start = cell->vertex(0); switch (dim) @@ -1148,6 +1186,8 @@ MappingCartesian::transform_real_to_unit_cell( const typename Triangulation::cell_iterator &cell, const Point & p) const { + Assert(is_cartesian(cell), ExcCellNotCartesian()); + if (dim != spacedim) Assert(false, ExcNotImplemented()); const Point &start = cell->vertex(0);