Jk/fix mapping

CMakeLists.txt

@@ -169,6 +169,7 @@ meshfields_add_exe(ElementJacobian2d test/testElementJacobian2d.cpp)
 meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
 meshfields_add_exe(OmegahElementTests test/testOmegahElement.cpp)
 meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
+meshfields_add_exe(PointMapping test/testPointMapping.cpp)
 
 if(MeshFields_USE_Cabana)
   meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
@@ -184,6 +185,7 @@ test_func(ElementJacobian1d ./ElementJacobian1d)
 test_func(ElementJacobian2d ./ElementJacobian2d)
 test_func(CountIntegrator ./CountIntegrator)
 test_func(OmegahElementTests ./OmegahElementTests)
+test_func(PointMapping ./PointMapping)
 if(MeshFields_USE_EXCEPTIONS)
   # exception caught - no error
   test_func(ExceptionTest ./ExceptionTest)

src/MeshField.hpp

@@ -57,7 +57,7 @@ createCoordinateField(const MeshField::MeshInfo &mesh_info,
 namespace MeshField {
 
 namespace Omegah {
-struct LinearTriangleToVertexField {
+template <size_t adjust = 0> struct LinearTriangleToVertexField {
   Omega_h::LOs triVerts;
   LinearTriangleToVertexField(Omega_h::Mesh &mesh)
       : triVerts(mesh.ask_elem_verts()) {
@@ -80,14 +80,23 @@ struct LinearTriangleToVertexField {
     const auto triDim = 2;
     const auto vtxDim = 0;
     const auto ignored = -1;
-    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
+    if (adjust == 1) {
+      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
+    } else {
+      const auto localVtxIdx =
+          (Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored) +
+           1) %
+          3;
+      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
+    }
   }
 };
-
 struct QuadraticTriangleToField {
   Omega_h::LOs triVerts;
   Omega_h::LOs triEdges;
@@ -153,10 +162,10 @@ template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
   if constexpr (ShapeOrder == 1) {
     struct result {
       MeshField::LinearTriangleShape shp;
-      LinearTriangleToVertexField map;
+      LinearTriangleToVertexField<1> map;
     };
     return result{MeshField::LinearTriangleShape(),
-                  LinearTriangleToVertexField(mesh)};
+                  LinearTriangleToVertexField<1>(mesh)};
   } else if constexpr (ShapeOrder == 2) {
     struct result {
       MeshField::QuadraticTriangleShape shp;

test/testOmegahElement.cpp

@@ -50,7 +50,7 @@ bool checkResult(Omega_h::Mesh &mesh, Result &result, CoordField coordField,
   const auto numPtsPerElem = testCase.NumPtsPerElem;
   MeshField::FieldElement fcoords(
       mesh.nfaces(), coordField, MeshField::LinearTriangleCoordinateShape(),
-      MeshField::Omegah::LinearTriangleToVertexField(mesh));
+      MeshField::Omegah::LinearTriangleToVertexField<1>(mesh));
   auto globalCoords =
       MeshField::evaluate(fcoords, testCase.coords, numPtsPerElem);
 

test/testPointMapping.cpp

@@ -0,0 +1,70 @@
+#include <Kokkos_Core.hpp>
+#include <MeshField.hpp>
+#include <Omega_h_build.hpp>
+#include <Omega_h_mesh.hpp>
+
+int main(int argc, char **argv) {
+  // setup
+  Kokkos::initialize(argc, argv);
+  auto lib = Omega_h::Library(&argc, &argv);
+  {
+    auto world = lib.world();
+    auto mesh =
+        Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false);
+    MeshField::OmegahMeshField<Kokkos::DefaultExecutionSpace, 2> mesh_field(
+        mesh);
+    auto shape_field = mesh_field.template CreateLagrangeField<double, 1, 1>();
+
+    // initialize vertices
+    auto dim = mesh.dim();
+    auto coords = mesh.coords();
+    auto f = KOKKOS_LAMBDA(double x, double y)->double {
+      return 2 * x + 3 * y;
+    };
+    Kokkos::parallel_for(
+        mesh.nverts(), KOKKOS_LAMBDA(int vtx) {
+          auto x = coords[vtx * dim];
+          auto y = coords[vtx * dim + 1];
+          shape_field(vtx, 0, 0, MeshField::Mesh_Topology::Vertex) = f(x, y);
+        });
+
+    Kokkos::View<double **> local_coords("", mesh.nelems() * 3, 3);
+    Kokkos::parallel_for(
+        "set", mesh.nelems() * 3, KOKKOS_LAMBDA(const int &i) {
+          local_coords(i, 0) = (i % 3 == 0);
+          local_coords(i, 1) = (i % 3 == 1);
+          local_coords(i, 2) = (i % 3 == 2);
+        });
+
+    auto eval_results =
+        mesh_field.triangleLocalPointEval(local_coords, 3, shape_field);
+
+    MeshField::Omegah::LinearTriangleToVertexField map(mesh);
+    int errors = 0;
+    Kokkos::parallel_reduce(
+        "test", mesh.nelems(),
+        KOKKOS_LAMBDA(const int &tri, int &errors) {
+          for (int node = 0; node < 3; ++node) {
+            auto mapping =
+                map(node, 0, tri, MeshField::Mesh_Topology::Triangle);
+            int vtx = mapping.entity;
+            auto x = coords[vtx * dim];
+            auto y = coords[vtx * dim + 1];
+            auto expected = f(x, y);
+            auto actual = eval_results(tri * 3 + node, 0);
+            if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) {
+              ++errors;
+              Kokkos::printf(
+                  "expected: %lf, actual: %lf, element: %d, node: %d\n",
+                  expected, actual, tri, node);
+            }
+          }
+        },
+        errors);
+    if (errors > 0) {
+      MeshField::fail("One or more mappings did not match\n");
+    }
+  }
+  Kokkos::finalize();
+  return 0;
+}