Question regarding FE<Dim,T>::reinit()

[nmnobre]

Hi,

Basically, for orientation-dependent element families such as Nédélec (first kind), even though shapes_need_reinit() returns true:

libmesh/src/fe/fe_nedelec_one.C

Lines 513 to 514 in 713688a

	template <> bool FE<2,NEDELEC_ONE>::shapes_need_reinit() const { return true; }
	template <> bool FE<3,NEDELEC_ONE>::shapes_need_reinit() const { return true; }

the following if statement:

libmesh/src/fe/fe.C

Lines 231 to 239 in 713688a

	if (this->shapes_need_reinit() && !cached_nodes_still_fit)
	{
	this->_fe_map->template init_reference_to_physical_map<Dim>
	(this->qrule->get_points(), elem);
	this->init_shape_functions (this->qrule->get_points(), elem);
	cached_nodes.resize(elem->n_nodes());
	for (auto n : elem->node_index_range())
	cached_nodes[n] = elem->point(n);
	}

will only execute if cached_nodes_still_fit is false. For sufficiently regular grids (simple compositions of quads or hexes), however, it will be true, resulting in no new computation. This implies that the computed orientation for the quads/hexes on the grid must also match that of the first element, is that why you use point coordinates instead of global ids to compute orientation?

Cheers,
Nuno (cc @karthichockalingam)

[roystgnr]

This implies that the computed orientation for the quads/hexes on the grid must also match that of the first element

Not ⇒, ⇐. If the orientations match, then cached_nodes_still_fit and we don't recompute. If they don't, we do. That's how it's supposed to work, anyway; you're right that we don't have nearly the level of test coverage we should here.

is that why you use point coordinates instead of global ids to compute orientation?

The main reason for avoiding global ids is that we like to be able to renumber those; otherwise doing transient adaptive mesh refinement+coarsening can leave you with a very sparse id set very quickly, and lots of codes (perhaps unwisely, but understandably) want to use ids as efficient global array indices.

On the other hand, the main reason to avoid using point coordinates is that doing so makes orientation-dependent finite elements unusable on problems with moving meshes. I could probably be talked into switching back (after rewriting our renumbering code to always be monotone; right now it does things like sorting by processor id on a DistributedMesh...).

[roystgnr]

Pinging @cticenhour since he might find the conversation interesting, either in a reassuring "yeah, we're using Nedelec on far more complicated meshes in the MOOSE EM module tests" way or in a dismaying "what do you mean we're not supposed to be using them on moving meshes" way...

[nmnobre]

Let me rephrase, the orientation needs to be computed using the lexicographic order on point coordinates and not the natural order on global ids, because otherwise, the reinit() code, which only checks the element spatial configuration to determine cached_nodes_still_fit, will incorrectly assume the orientation (as given by the global ids) is also preserved.

[roystgnr]

I sure let this get buried, didn't I? I think we came to an understanding, despite me misunderstanding your question at first, so I'll close the issue unless there's still something confusing.

otherwise, the reinit() code, which only checks the element spatial configuration to determine cached_nodes_still_fit, will incorrectly assume the orientation (as given by the global ids) is also preserved.

This is indeed the case. Which means that as soon as we add an option to determine orientation via unique_id() rather than geometry, we're also going to need to ... probably disable caching entirely, when an element needs that option? "This element is just a translation of the last one you computed" is an extremely common case, but "this element also happens to be numbered with the exact same ordering" is a 1/N! case, not even worth testing for. On the bright side, the displaced-mesh cases where we can't trust geometrically determined orderings to stay unchanged are also the cases where we're perturbing every element and can't cache much anyway.