Kokkos support for Lagrange and Monomial elements

[rochi00]

1. Shape functions — nativeShape, nativeGradShape for LAGRANGE (12 types) + MONOMIAL (orders 0-5)
2. Type dispatch — FEShapeKey, classFromTopology, nDofs, templated nativeMapShape
3. Quadrature — GaussQuadrature device-callable for EDGE, QUAD, HEX, TRI, TET + fillQuadrature host helper
4. Physical map — physicalPoint, jacobian, physicalPointAndJacobian, faceJacobian, faceJxW, faceNormal (templated + runtime versions)
5. Device-safe asserts — printf + Kokkos::abort() for debug builds on GPU

[moosebuild]

Job Coverage, step Generate coverage on 639e20b wanted to post the following:

Coverage

	325622	#4441 639e20
	Total	Total	+/-	New
Rate	65.47%	65.47%	+0.00%	100.00%
Hits	78228	78230	+2	6
Misses	41258	41256	-2	0

Diff coverage report

Full coverage report

This comment will be updated on new commits.

[roystgnr]

In practice I think literally every compiler we care about supports this (though Wiki lists one I've never heard of that doesn't!), but since it never made it into a C++ standard we've always used include guards instead; let's do that here too.

[roystgnr]

Mostly when we throw an exception it's because we can't sanely go on; even if we don't support exceptions on device I don't think silently going on is the right way to handle that. When we don't support exceptions on the host we print the error message and abort; should we not do that here too?

[roystgnr]

Concept is good, but I don't like the name. FEFamily really is FE Method specific; an Elem is an Elem even if we're doing Finite Volume or pure geometry/meshing work (or in the future FDM stencils or whatever) instead. Let's just call this ElemClass, unless someone has an idea for a better name?

[roystgnr]

No it doesn't. p=1 Lagrange on a Quad9 makes perfect sense and is commonly used for the pressure variable in incompressible flow simulation.

I get that we're not supporting everything right out of the gate, but we need to make sure that we're not digging ourselves into any holes in our documentation or (more importantly, since it'll be much harder to change later) our APIs. Is there something in the design that would already prevent us from returning evaluations of a p=1 Lagrange on a higher-order geometric element?

[roystgnr]

I'd like to use snake_case rather than camelCase, like we do everywhere else in libMesh code.

We should also try to get this to resemble existing code in less superficial ways if possible. We've got static const unsigned int Elem::type_to_dim_map[] already, for instance; could we either add an Elem::class_to_dim_map[] too, or if elem.h would cause Kokkos fits, would it make sense for us to still use the same array-based design to make things easier to merge later?

[roystgnr]

They don't, because that's kind of horrifying for type safety. Even if they did we'd probably decide that the proper projection of s from ℝ→ℝ³ was (s,0,0), not (s,s,s). That's not even of magnitude s!

Where do we use these? If the answer is "nowhere" let's delete them; if not let's give them new names so I can grep for those names and stare at the hits suspiciously.

[roystgnr]

Wait, what?

There are reference coordinates on a point element; they're all just (0,0,0).

But (x,y,z) is what point(0) will give you, and that is not (xi,eta,zeta). point(0) on a NodeElem will be a point in physical space.

[roystgnr]

Okay, it's midnight, and I'm getting tired of tearing my hair out, and I think I've already written enough to keep you going for a while, so I'm going to post this review half done, with one final issue:

Is this whole block not just wrong? It gives the location in physical space of face_qpt in side-reference-space, but the docs here say we're returning a point in parent-reference space. There has to be an inverse_map() afterward to get back to that.

Please tell me I'm wrong, and I just lose too many IQ points by midnight? We clearly need way more test coverage if I'm somehow our first line of defense for a bug this big. I know we don't have the environments set up to test in PR yet, but ideally you're running make check yourself and at some point in the new tests we're doing a side integral?

[roystgnr]

No. Definitely not.

[roystgnr]

Actually, one more comment.

There's an interesting story about the rock band Van Halen:

When they were on tour, they had a fairly complicated and safety-critical set of contract stipulations that venues had to meet: heavy platforms and heavy equipment, serious electrical requirements, pyrotechnics, etc. And they had a little contract clause that demanded a bowl of M&Ms for them backstage, with all the brown M&Ms removed.

If someone challenged that clause, they waived it. They didn't care about not eating brown M&Ms. They just cared that the big long critical pile of contract text had actually been respected, entirely read by the people with the primary responsibility of reading it and making sure everything was set up correctly. If they hadn't been reading carefully enough to notice the brown-M&M-hating-divas bit, could they really be trusted to be getting the subtle, tricky mechanical and electrical and thermal stuff done and inspected safely? No; the band roadies would essentially have to reinspect everything to see which parts of it they'd have to redo.

It was not a good idea to be a venue that made the hard rock band and their annoyed roadies unable to trust any of your work.

[rochi00]

Sorry, I mostly paid attention to the tests that were written, and once they all passed assumed the internal code was correct. That was clearly not a good idea. I will review everything now. The tests clearly didn't have enough coverage

[rochi00]

The claude generated comments had so many issues, it is terrifying how much it overclaims