FE pre-requests for error estimation

[nmnobre]

Hi @jwpeterson,

Regarding the warning messages you're seeing in the newest vector examples:

As discussed over slack with @roystgnr, we may have a missing FE pre-request somewhere which is triggering the following "not currently being computed" warnings for me. Something to look into for a separate PR...

WARNING: Shape function gradients for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 88, compiled Oct 31 2023 at 13:52:00 ***
WARNING: Shape function curls for HDiv elements are not currently being computed!./include/libmesh/hdiv_fe_transformation.h, line 121, compiled Oct 31 2023 at 13:52:00 ***

Originally posted by @jwpeterson in #3686 (review)

The culprit is the following piece of code used to evaluate the computed solution for error estimation purposes:

libmesh/src/error_estimation/exact_solution.C

Lines 616 to 653 in b2248b5

	// The value of the shape functions at the quadrature points
	// i.e. phi(i) = phi_values[i][qp]
	const std::vector<std::vector<OutputShape>> & phi_values = fe->get_phi();
	// The value of the shape function gradients at the quadrature points
	const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputGradient>> &
	dphi_values = fe->get_dphi();
	// The value of the shape function curls at the quadrature points
	// Only computed for vector-valued elements
	const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputShape>> * curl_values = nullptr;
	// The value of the shape function divergences at the quadrature points
	// Only computed for vector-valued elements
	const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputDivergence>> * div_values = nullptr;
	if (field_type == TYPE_VECTOR)
	{
	curl_values = &fe->get_curl_phi();
	div_values = &fe->get_div_phi();
	}
	#ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
	// The value of the shape function second derivatives at the quadrature points
	// Not computed for vector-valued elements
	const std::vector<std::vector<typename FEGenericBase<OutputShape>::OutputTensor>> *
	d2phi_values = nullptr;
	if (field_type != TYPE_VECTOR)
	d2phi_values = &fe->get_d2phi();
	#endif
	// The XYZ locations (in physical space) of the quadrature points
	const std::vector<Point> & q_point = fe->get_xyz();
	// reinitialize the element-specific data
	// for the current element
	fe->reinit (elem);

We could, perhaps, add a couple of guards, or maybe just implement the missing maps, even if they aren't really used anywhere xD

[lindsayad]

nice sleuthing!

[jwpeterson]

OK, so it's nothing to do with missing pre-requests, it's that the call to fe->get_dphi() calls HDivFETransformation::map_dphi() which is unimplemented?

even if they aren't really used anywhere

I'm not sure what you mean, dphi_values is definitely used in the ExactSolution function that you highlighted?

[nmnobre]

it's that the call to fe->get_dphi() calls HDivFETransformation::map_dphi() which is unimplemented

Indeed.

I'm not sure what you mean, dphi_values is definitely used in the ExactSolution function that you highlighted?

In the sense that, for instance for ex6/ex7, you ask for gradients but don't compute errors w/ the H1 norm, you ask for curls, but don't compute errors w/ the H(curl) norm. That's why the warnings can be safely ignored in those cases, they just don't look very nice.

So what I meant was, we can implement the maps, and they'll be used at runtime, as you rightly said.
But unless the user asks for, e.g. H1/H(curl) errors for H(div) elements, they just won't affect the semantics of the program, it'll be dead code (as it currently is, it's just all zeros).