shape functions: order of parametric coords

[cwsmith]

meshFields/src/MeshField_Shape.hpp

Lines 75 to 77 in cb7fc5e

	return {1 - xi[0] - xi[1],
	xi[0],
	xi[1]};

Why are we using the 3rd component first?

[cwsmith]

This follows what was done in pumi:

    class Triangle : public EntityShape
    {    
      public:
        void getValues(Mesh*, MeshEntity*,
            Vector3 const& xi, NewArray<double>& values) const
        {    
          values.allocate(3);
          values[0] = 1-xi[0]-xi[1];
          values[1] = xi[0];
          values[2] = xi[1];

from here: https://github.com/SCOREC/core/blob/18eda751efad04d095e9da318a8244e9e8d4dafd/apf/apfShape.cc#L140-L149

From what I understand, any order is fine as long as the user understands the convention and the code is internally consistent. I'm not sure why it was written this way in pumi.

I'm not opposed to changing this if that would be more intuitive. @jacobmerson

[jacobmerson]

I do think that the other way seems more intuitive because then the parametric coordinates are the same as the barycentric coordinates for the linear case. But, we will need to make sure we are careful that the node ordering is matching in the integration rules as well.

We may want to confirm with Mark if he remembers a specific reason why this ordering was chosen. It's not obvious to me, but that doesn't mean it wasn't done for a good reason.

[cwsmith]

Mark didn't know offhand why the ordering was that way.

[cwsmith]

I've found a couple references that use the same convention as pumi:

• https://users.metu.edu.tr/csert/me582/ME582%20Ch%2003.pdf
• https://pkel015.connect.amazon.auckland.ac.nz/SolidMechanicsBooks/FEM/Two_Dimensional/07_Triangles.pdf

and MFEM does the same (https://docs.mfem.org/4.8/fe__fixed__order_8cpp_source.html#l00074). Given that, I don't think we should change the convention, but need to document it.

[joyal-mathew]

I believe there is an inconsistency between this and the field accessors. LinearTriangleToVertexField returns vertices in a different order than here so when evaluating the field, the wrong coefficients are applied to the nodal values.

Right now, I'm reordering the coordinates before I pass them into meshFields to get the correct result.

[cwsmith]

If you can you provide a small reproducer that would be helpful.

Note to self - look into the indexing here:

meshFields/src/MeshField.hpp

Lines 83 to 87 in 29104f8

	const auto localVtxIdx =
	Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored);
	const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim);
	const MeshField::LO vtx = triVerts[(tri * triToVtxDegree) + localVtxIdx];
	return {0, triCompIdx, vtx, MeshField::Vertex}; // node, comp, ent, topo

[joyal-mathew]

Here's an example of what I mean:
https://gist.github.com/joyal-mathew/d5a20a6a48e36bd9e214ec8df4831094

I explain more in the comments, but essentially, I would expect parametric coordinate 0 and node 0 (when using LinearTriangleToVertexField) to refer to the same vertex of a given triangle. But this is not the case, which causes an issue with evaluation.

[cwsmith]

@joyal-mathew Great. Thank you.

[cwsmith]

@joyal-mathew Thanks again for reporting this.
#64 was just merged that resolves the ordering issue you hit. Please take a look at this test for example usage using the omegah canonical ordering apis and let @Joshua-Kloepfer and I know if there are questions/concerns/etc.

meshFields/test/testPointMapping.cpp

Lines 64 to 82 in 05fc74d

	for (int node = 0; node < 3; ++node) {
	const auto triDim = 2;
	const auto vtxDim = 0;
	const auto ignored = -1;
	const auto localVtxIdx =
	Omega_h::simplex_down_template(triDim, vtxDim, node, ignored);
	const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim);
	int vtx = triVerts[(tri * triToVtxDegree) + localVtxIdx];
	auto x = coords[vtx * dim];
	auto y = coords[vtx * dim + 1];
	auto expected = f(x, y);
	auto actual = eval_results(tri * numNodesPerElem + node, 0);
	if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) {
	++errors;
	Kokkos::printf(
	"expected: %lf, actual: %lf, element: %d, node(vtx): %d\n",
	expected, actual, tri, node);
	}
	}