From 3f901b374d4118fdf38209cc1c044f398782cc4b Mon Sep 17 00:00:00 2001 From: Kloepfer Date: Thu, 17 Jul 2025 17:17:23 -0400 Subject: [PATCH 1/6] Working on test for mappings. Co-authored-by: Joyal Mathew --- CMakeLists.txt | 2 + src/MeshField.hpp | 27 +++++++++----- test/testOmegahElement.cpp | 2 +- test/testPointMapping.cpp | 75 ++++++++++++++++++++++++++++++++++++++ 4 files changed, 96 insertions(+), 10 deletions(-) create mode 100644 test/testPointMapping.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d3b6ef..b78082b 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 0be91e8..5466cf2 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -57,7 +57,7 @@ createCoordinateField(const MeshField::MeshInfo &mesh_info, namespace MeshField { namespace Omegah { -struct LinearTriangleToVertexField { +template 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 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; diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index 8e055ed..b334918 100644 --- a/test/testOmegahElement.cpp +++ b/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); diff --git a/test/testPointMapping.cpp b/test/testPointMapping.cpp new file mode 100644 index 0000000..eaca8eb --- /dev/null +++ b/test/testPointMapping.cpp @@ -0,0 +1,75 @@ +#include +#include +#include +#include + +// Uncommenting will reorder the parametric evaluation coordinates which causes +// the test to pass. +// #define REORDER_COORDS + +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 mesh_field( + mesh); + auto shape_field = mesh_field.template CreateLagrangeField(); + + // initialize vertices + auto dim = mesh.dim(); + auto coords = mesh.coords(); + auto f = KOKKOS_LAMBDA(double x, double y)->double { + return 2 * x + 3 * y; + }; + // auto f = [](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 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); + 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; +} From d22d60406d113301d4ca6152dd6d5de3b8786b94 Mon Sep 17 00:00:00 2001 From: Kloepfer Date: Fri, 18 Jul 2025 12:17:30 -0400 Subject: [PATCH 2/6] Fixed typo Co-authored-by: Joyal Mathew --- test/testPointMapping.cpp | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/test/testPointMapping.cpp b/test/testPointMapping.cpp index eaca8eb..ed62652 100644 --- a/test/testPointMapping.cpp +++ b/test/testPointMapping.cpp @@ -3,10 +3,6 @@ #include #include -// Uncommenting will reorder the parametric evaluation coordinates which causes -// the test to pass. -// #define REORDER_COORDS - int main(int argc, char **argv) { // setup Kokkos::initialize(argc, argv); @@ -25,7 +21,6 @@ int main(int argc, char **argv) { auto f = KOKKOS_LAMBDA(double x, double y)->double { return 2 * x + 3 * y; }; - // auto f = [](double x, double y) -> double { return 2 * x + 3 * y; }; Kokkos::parallel_for( mesh.nverts(), KOKKOS_LAMBDA(int vtx) { auto x = coords[vtx * dim]; @@ -56,7 +51,7 @@ int main(int argc, char **argv) { auto x = coords[vtx * dim]; auto y = coords[vtx * dim + 1]; auto expected = f(x, y); - auto actual = eval_results(tri * 3, node); + auto actual = eval_results(tri * 3 + node, 0); if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) { ++errors; Kokkos::printf( From 5c0313e33f38e6400f36376729041e5c7cf06f40 Mon Sep 17 00:00:00 2001 From: Kloepfer Date: Tue, 22 Jul 2025 11:46:17 -0400 Subject: [PATCH 3/6] Adjusted tests --- src/MeshField.hpp | 22 ++++++------------ test/testOmegahElement.cpp | 2 +- test/testPointMapping.cpp | 46 ++++++++++++++++++++++---------------- 3 files changed, 35 insertions(+), 35 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 5466cf2..5337d70 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -57,7 +57,7 @@ createCoordinateField(const MeshField::MeshInfo &mesh_info, namespace MeshField { namespace Omegah { -template struct LinearTriangleToVertexField { +struct LinearTriangleToVertexField { Omega_h::LOs triVerts; LinearTriangleToVertexField(Omega_h::Mesh &mesh) : triVerts(mesh.ask_elem_verts()) { @@ -80,21 +80,13 @@ template struct LinearTriangleToVertexField { const auto triDim = 2; const auto vtxDim = 0; const auto ignored = -1; - 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 = + const auto localVtxIdx = (Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored) + - 1) % + 2) % 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 { @@ -139,7 +131,7 @@ struct QuadraticTriangleToField { const auto vtxDim = 0; const auto ignored = -1; const auto localVtxIdx = - Omega_h::simplex_down_template(triDim, vtxDim, dofHolderIdx, ignored); + (Omega_h::simplex_down_template(triDim, vtxDim, dofHolderIdx, ignored) + 2) % 3; const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim); osh_ent = triVerts[(tri * triToVtxDegree) + localVtxIdx]; } else if (dofHolderTopo == MeshField::Edge) { @@ -149,7 +141,7 @@ struct QuadraticTriangleToField { // passing dofHolderIdx as Omega_h_simplex.hpp does not provide // a function that maps a triangle and edge index to a 'canonical' edge // index. This may need to be revisited... - osh_ent = triEdges[(tri * triToEdgeDegree) + dofHolderIdx]; + osh_ent = triEdges[(tri * triToEdgeDegree) + (dofHolderIdx + 2) % 3]; } else { assert(false); } @@ -162,10 +154,10 @@ template auto getTriangleElement(Omega_h::Mesh &mesh) { if constexpr (ShapeOrder == 1) { struct result { MeshField::LinearTriangleShape shp; - LinearTriangleToVertexField<1> map; + LinearTriangleToVertexField map; }; return result{MeshField::LinearTriangleShape(), - LinearTriangleToVertexField<1>(mesh)}; + LinearTriangleToVertexField(mesh)}; } else if constexpr (ShapeOrder == 2) { struct result { MeshField::QuadraticTriangleShape shp; diff --git a/test/testOmegahElement.cpp b/test/testOmegahElement.cpp index b334918..8e055ed 100644 --- a/test/testOmegahElement.cpp +++ b/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<1>(mesh)); + MeshField::Omegah::LinearTriangleToVertexField(mesh)); auto globalCoords = MeshField::evaluate(fcoords, testCase.coords, numPtsPerElem); diff --git a/test/testPointMapping.cpp b/test/testPointMapping.cpp index ed62652..9d62cc4 100644 --- a/test/testPointMapping.cpp +++ b/test/testPointMapping.cpp @@ -2,18 +2,11 @@ #include #include #include - -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); +template +void runTest(Omega_h::Mesh mesh) { MeshField::OmegahMeshField mesh_field( mesh); - auto shape_field = mesh_field.template CreateLagrangeField(); + auto shape_field = mesh_field.template CreateLagrangeField(); // initialize vertices auto dim = mesh.dim(); @@ -27,8 +20,8 @@ int main(int argc, char **argv) { auto y = coords[vtx * dim + 1]; shape_field(vtx, 0, 0, MeshField::Mesh_Topology::Vertex) = f(x, y); }); - - Kokkos::View local_coords("", mesh.nelems() * 3, 3); + const auto numPtsPerElem = 3; + Kokkos::View local_coords("", mesh.nelems() * numPtsPerElem, 3); Kokkos::parallel_for( "set", mesh.nelems() * 3, KOKKOS_LAMBDA(const int &i) { local_coords(i, 0) = (i % 3 == 0); @@ -37,21 +30,24 @@ int main(int argc, char **argv) { }); auto eval_results = - mesh_field.triangleLocalPointEval(local_coords, 3, shape_field); + mesh_field.triangleLocalPointEval(local_coords, numPtsPerElem, shape_field); - MeshField::Omegah::LinearTriangleToVertexField map(mesh); int errors = 0; + const auto triVerts = mesh.ask_elem_verts(); 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; + for (int node = 0; node < numPtsPerElem; ++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 * 3 + node, 0); + auto actual = eval_results(tri * numPtsPerElem + node, 0); if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) { ++errors; Kokkos::printf( @@ -64,6 +60,18 @@ int main(int argc, char **argv) { if (errors > 0) { MeshField::fail("One or more mappings did not match\n"); } + +} +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); + runTest<1>(mesh); + runTest<2>(mesh); } Kokkos::finalize(); return 0; From c5c486f9468e272c18f75f616db246a3d66f4f17 Mon Sep 17 00:00:00 2001 From: Kloepfer Date: Tue, 22 Jul 2025 11:50:29 -0400 Subject: [PATCH 4/6] formatting --- src/MeshField.hpp | 20 +++---- test/testPointMapping.cpp | 107 ++++++++++++++++++-------------------- 2 files changed, 63 insertions(+), 64 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 5337d70..dd94241 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -80,13 +80,13 @@ 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) + - 2) % - 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 + const auto localVtxIdx = + (Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored) + + 2) % + 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 { @@ -130,8 +130,10 @@ struct QuadraticTriangleToField { const auto triDim = 2; const auto vtxDim = 0; const auto ignored = -1; - const auto localVtxIdx = - (Omega_h::simplex_down_template(triDim, vtxDim, dofHolderIdx, ignored) + 2) % 3; + const auto localVtxIdx = (Omega_h::simplex_down_template( + triDim, vtxDim, dofHolderIdx, ignored) + + 2) % + 3; const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim); osh_ent = triVerts[(tri * triToVtxDegree) + localVtxIdx]; } else if (dofHolderTopo == MeshField::Edge) { diff --git a/test/testPointMapping.cpp b/test/testPointMapping.cpp index 9d62cc4..7b37f55 100644 --- a/test/testPointMapping.cpp +++ b/test/testPointMapping.cpp @@ -2,65 +2,62 @@ #include #include #include -template -void runTest(Omega_h::Mesh mesh) { - MeshField::OmegahMeshField mesh_field( - mesh); - auto shape_field = mesh_field.template CreateLagrangeField(); +template void runTest(Omega_h::Mesh mesh) { + MeshField::OmegahMeshField mesh_field(mesh); + auto shape_field = + mesh_field.template CreateLagrangeField(); - // 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); - }); - const auto numPtsPerElem = 3; - Kokkos::View local_coords("", mesh.nelems() * numPtsPerElem, 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); - }); + // 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); + }); + const auto numPtsPerElem = 3; + Kokkos::View local_coords("", mesh.nelems() * numPtsPerElem, 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, numPtsPerElem, shape_field); + auto eval_results = mesh_field.triangleLocalPointEval( + local_coords, numPtsPerElem, shape_field); - int errors = 0; - const auto triVerts = mesh.ask_elem_verts(); - Kokkos::parallel_reduce( - "test", mesh.nelems(), - KOKKOS_LAMBDA(const int &tri, int &errors) { - for (int node = 0; node < numPtsPerElem; ++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 * numPtsPerElem + 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); - } + int errors = 0; + const auto triVerts = mesh.ask_elem_verts(); + Kokkos::parallel_reduce( + "test", mesh.nelems(), + KOKKOS_LAMBDA(const int &tri, int &errors) { + for (int node = 0; node < numPtsPerElem; ++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 * numPtsPerElem + 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"); - } - + } + }, + errors); + if (errors > 0) { + MeshField::fail("One or more mappings did not match\n"); + } } int main(int argc, char **argv) { // setup From 204ec428ad542cc7de7a5acc5d911816a6c06d9d Mon Sep 17 00:00:00 2001 From: Kloepfer Date: Fri, 25 Jul 2025 12:21:05 -0400 Subject: [PATCH 5/6] added edge check for quadratic case --- test/testPointMapping.cpp | 63 ++++++++++++++++++++++++++++++++++----- 1 file changed, 55 insertions(+), 8 deletions(-) diff --git a/test/testPointMapping.cpp b/test/testPointMapping.cpp index 7b37f55..728d15b 100644 --- a/test/testPointMapping.cpp +++ b/test/testPointMapping.cpp @@ -11,19 +11,45 @@ template void runTest(Omega_h::Mesh mesh) { auto dim = mesh.dim(); auto coords = mesh.coords(); auto f = KOKKOS_LAMBDA(double x, double y)->double { return 2 * x + 3 * y; }; + auto edge2vtx = mesh.ask_down(1, 0).ab2b; + auto edgeMap = mesh.ask_down(dim, 1).ab2b; 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); }); - const auto numPtsPerElem = 3; + if constexpr (order == 2) { + Kokkos::parallel_for( + mesh.nedges(), KOKKOS_LAMBDA(int edge) { + const auto left = edge2vtx[edge * 2]; + const auto right = edge2vtx[edge * 2 + 1]; + const auto x = (coords[left * dim] + coords[right * dim]) / 2.0; + const auto y = + (coords[left * dim + 1] + coords[right * dim + 1]) / 2.0; + shape_field(edge, 0, 0, MeshField::Mesh_Topology::Edge) = f(x, y); + }); + } + const auto numPtsPerElem = order == 2 ? 6 : 3; Kokkos::View local_coords("", mesh.nelems() * numPtsPerElem, 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); + "set", mesh.nelems() * numPtsPerElem, KOKKOS_LAMBDA(const int &i) { + const auto val = i % numPtsPerElem; + local_coords(i, 0) = (val == 0); + local_coords(i, 1) = (val == 1); + local_coords(i, 2) = (val == 2); + if constexpr (order == 2) { + if (val == 3) { + local_coords(i, 0) = 1 / 2.0; + local_coords(i, 1) = 1 / 2.0; + } else if (val == 4) { + local_coords(i, 1) = 1 / 2.0; + local_coords(i, 2) = 1 / 2.0; + } else if (val == 5) { + local_coords(i, 0) = 1 / 2.0; + local_coords(i, 2) = 1 / 2.0; + } + } }); auto eval_results = mesh_field.triangleLocalPointEval( @@ -31,10 +57,11 @@ template void runTest(Omega_h::Mesh mesh) { int errors = 0; const auto triVerts = mesh.ask_elem_verts(); + const auto coordField = mesh_field.getCoordField(); Kokkos::parallel_reduce( "test", mesh.nelems(), KOKKOS_LAMBDA(const int &tri, int &errors) { - for (int node = 0; node < numPtsPerElem; ++node) { + for (int node = 0; node < 3; ++node) { const auto triDim = 2; const auto vtxDim = 0; const auto ignored = -1; @@ -49,8 +76,28 @@ template void runTest(Omega_h::Mesh mesh) { if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) { ++errors; Kokkos::printf( - "expected: %lf, actual: %lf, element: %d, node: %d\n", expected, - actual, tri, node); + "expected: %lf, actual: %lf, element: %d, node(vtx): %d\n", + expected, actual, tri, node); + } + } + for (int node = 3; node < numPtsPerElem; ++node) { + const auto triDim = 2; + const auto edgeDim = 1; + const auto triToEdgeDegree = Omega_h::simplex_degree(triDim, edgeDim); + const MeshField::LO triNode2DofHolder[3] = {0, 1, 2}; + int edge = + edgeMap[(tri * triToEdgeDegree) + triNode2DofHolder[node - 3]]; + auto left = edge2vtx[edge * 2]; + auto right = edge2vtx[edge * 2 + 1]; + auto x = (coords[left * dim] + coords[right * dim]) / 2.0; + auto y = (coords[left * dim + 1] + coords[right * dim + 1]) / 2.0; + auto expected = f(x, y); + auto actual = eval_results(tri * numPtsPerElem + node, 0); + if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) { + ++errors; + Kokkos::printf( + "expected: %lf, actual: %lf, element: %d, node(edge): %d\n", + expected, actual, tri, node); } } }, From d00c854bad5277a39c297b25192a876a11097031 Mon Sep 17 00:00:00 2001 From: Kloepfer Date: Fri, 25 Jul 2025 14:03:30 -0400 Subject: [PATCH 6/6] fixes --- test/testPointMapping.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/test/testPointMapping.cpp b/test/testPointMapping.cpp index 728d15b..e71c8e0 100644 --- a/test/testPointMapping.cpp +++ b/test/testPointMapping.cpp @@ -19,7 +19,7 @@ template void runTest(Omega_h::Mesh mesh) { auto y = coords[vtx * dim + 1]; shape_field(vtx, 0, 0, MeshField::Mesh_Topology::Vertex) = f(x, y); }); - if constexpr (order == 2) { + if (order == 2) { Kokkos::parallel_for( mesh.nedges(), KOKKOS_LAMBDA(int edge) { const auto left = edge2vtx[edge * 2]; @@ -30,11 +30,11 @@ template void runTest(Omega_h::Mesh mesh) { shape_field(edge, 0, 0, MeshField::Mesh_Topology::Edge) = f(x, y); }); } - const auto numPtsPerElem = order == 2 ? 6 : 3; - Kokkos::View local_coords("", mesh.nelems() * numPtsPerElem, 3); + const auto numNodesPerElem = order == 2 ? 6 : 3; + Kokkos::View local_coords("", mesh.nelems() * numNodesPerElem, 3); Kokkos::parallel_for( - "set", mesh.nelems() * numPtsPerElem, KOKKOS_LAMBDA(const int &i) { - const auto val = i % numPtsPerElem; + "set", mesh.nelems() * numNodesPerElem, KOKKOS_LAMBDA(const int &i) { + const auto val = i % numNodesPerElem; local_coords(i, 0) = (val == 0); local_coords(i, 1) = (val == 1); local_coords(i, 2) = (val == 2); @@ -53,7 +53,7 @@ template void runTest(Omega_h::Mesh mesh) { }); auto eval_results = mesh_field.triangleLocalPointEval( - local_coords, numPtsPerElem, shape_field); + local_coords, numNodesPerElem, shape_field); int errors = 0; const auto triVerts = mesh.ask_elem_verts(); @@ -72,7 +72,7 @@ template void runTest(Omega_h::Mesh mesh) { auto x = coords[vtx * dim]; auto y = coords[vtx * dim + 1]; auto expected = f(x, y); - auto actual = eval_results(tri * numPtsPerElem + node, 0); + auto actual = eval_results(tri * numNodesPerElem + node, 0); if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) { ++errors; Kokkos::printf( @@ -80,7 +80,7 @@ template void runTest(Omega_h::Mesh mesh) { expected, actual, tri, node); } } - for (int node = 3; node < numPtsPerElem; ++node) { + for (int node = 3; node < numNodesPerElem; ++node) { const auto triDim = 2; const auto edgeDim = 1; const auto triToEdgeDegree = Omega_h::simplex_degree(triDim, edgeDim); @@ -92,7 +92,7 @@ template void runTest(Omega_h::Mesh mesh) { auto x = (coords[left * dim] + coords[right * dim]) / 2.0; auto y = (coords[left * dim + 1] + coords[right * dim + 1]) / 2.0; auto expected = f(x, y); - auto actual = eval_results(tri * numPtsPerElem + node, 0); + auto actual = eval_results(tri * numNodesPerElem + node, 0); if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) { ++errors; Kokkos::printf(