MatrixFree: check that JxW is not negative

[peterrum]

In a follow-up PR, I would like to extend it like

dealii/source/fe/mapping_q.cc

Lines 932 to 936 in 9042b92

	Assert(det >
	1e-12 * Utilities::fixed_power<dim>(
	cell->diameter() / std::sqrt(double(dim))),
	(typename Mapping<dim, spacedim>::ExcDistortedMappedCell(
	cell->center(), det, point)));

But for that I access to the cells.

[peterrum]

@nfehn Could you check if this assert is triggered for your mesh?

[nfehn]

I think this is a very helpful check.

@nfehn Could you check if this assert is triggered for your mesh?

Yes, it is.

[nfehn]

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.

[kronbichler]

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.

[peterrum]

@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.

[kronbichler]

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.

[nfehn]

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.

[bangerth]

/rebuild

[peterrum]

@kronbichler Do you have an idea why matrix_free/matrix_vector_hessians_general_faces.debug is triggering the assert? matrix_free/matrix_vector_hessians_general_cells.debug is not. The only difference is that one sets mapping_update_flags_inner_faces and the other mapping_update_flags.