Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

MatrixFree: check that JxW is not negative #16317

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
+35 −4
Open

MatrixFree: check that JxW is not negative #16317

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9a3824e
MatrixFree: check that JxW is not negative
peterrum Dec 3, 2023
9356cc7
Extend assert message
peterrum Dec 3, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
39 changes: 35 additions & 4 deletions include/deal.II/matrix_free/mapping_info.templates.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1152,10 +1152,12 @@ namespace internal
typename VectorizedDouble>
void
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When looking at the code, I realized that the documentation for the other function below (dealing with faces), see comment https://github.com/dealii/dealii/pull/16317/files#r1413627125, seems to be the same as here for the cells? This could be something one wants to change in the future. Of course this is somewhat unrelated to the present PR.

mapping_q_compute_range(
const unsigned int begin_cell,
const unsigned int end_cell,
const std::vector<GeometryType> &cell_type,
const std::vector<bool> &process_cell,
const unsigned int begin_cell,
const unsigned int end_cell,
const dealii::Triangulation<dim> &tria,
const std::vector<std::pair<unsigned int, unsigned int>> &cell_array,
const std::vector<GeometryType> &cell_type,
const std::vector<bool> &process_cell,
const UpdateFlags update_flags_cells,
const AlignedVector<double> &plain_quadrature_points,
const ShapeInfo<VectorizedDouble> &shape_info,
Expand Down Expand Up @@ -1240,6 +1242,28 @@ namespace internal
jac[d][e] = 0.;

const VectorizedDouble jac_det = determinant(jac);

#ifdef DEBUG
for (unsigned int v = 0; v < n_lanes_d; ++v)
{
const typename Triangulation<dim>::cell_iterator
cell_iterator(
&tria,
cell_array[cell * n_lanes + vv + v].first,
cell_array[cell * n_lanes + vv + v].second);

Assert(jac_det[v] >
1e-12 * Utilities::fixed_power<dim>(
cell_iterator->diameter() /
std::sqrt(double(dim))),
(typename Mapping<dim>::ExcDistortedMappedCell(
cell_iterator->center(), jac_det[v], q)));
}
#else
(void)tria;
(void)cell_array;
#endif

const Tensor<2, dim, VectorizedDouble> inv_jac =
transpose(invert(jac));

Expand Down Expand Up @@ -2309,6 +2333,11 @@ namespace internal
1. :
my_data.descriptor[0].quadrature.weight(q));

#ifdef DEBUG
for (unsigned int v = 0; v < n_lanes_d; ++v)
Assert(JxW[v] > 0.0, ExcInternalError());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I wondered whether this exc message is clear, i.e. points to an invalid Jacobian determinant on faces? One might argue that JxW is basic FE language, so this is probably a matter of taste.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think this message should ever trigger, as opposed to the one in line 1254 above, because JxW is computed as

  JxW = boundary_form.norm() * my_data.descriptor[0].quadrature.weight(q);

so here we have a norm that cannot be negative and a quadrature weight that should be non-negative. So here I would guess it is ok.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@kronbichler Thank you for the explanation! But the value could be 0, couldn't it? Furthermore, the check for the cell not only checks whether the value is positive but also that it is not too small. Unfortunately, I could not find such a check for faces in MappingQ (maybe I missed it).

Should we remove this assert? The chance is quite high that assert for the cell is triggered.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think there is a check for the faces in terms of "sufficient" positivity. I would personally choose something like the other assertion you introduced above, with the difference to use Utilities::fixed_power<dim-1>(edge_length), i.e., one power less. It can trigger for really anisotropic faces (aspect ratio larger than 1e5), but I think that is ok.

#endif

store_vectorized_array(JxW,
vv,
my_data.JxW_values[offset + q]);
Expand Down Expand Up @@ -2793,6 +2822,8 @@ namespace internal
VectorizedDouble>(
begin,
end,
tria,
cell_array,
cell_type,
process_cell,
update_flags_cells,
Expand Down