Skip to content

Commit 05fc74d

Browse files
cwsmithcwsmith
authored
Merge pull request #64 from SCOREC/jk/fixMapping
Jk/fix mapping
2 parents b9319c5 + d00c854 commit 05fc74d

File tree

3 files changed

+132
-5
lines changed

3 files changed

+132
-5
lines changed

‎CMakeLists.txt‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -169,6 +169,7 @@ meshfields_add_exe(ElementJacobian2d test/testElementJacobian2d.cpp)
169169
meshfields_add_exe(CountIntegrator test/testCountIntegrator.cpp)
170170
meshfields_add_exe(OmegahElementTests test/testOmegahElement.cpp)
171171
meshfields_add_exe(ExceptionTest test/testExceptions.cpp)
172+
meshfields_add_exe(PointMapping test/testPointMapping.cpp)
172173

173174
if(MeshFields_USE_Cabana)
174175
meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp)
@@ -184,6 +185,7 @@ test_func(ElementJacobian1d ./ElementJacobian1d)
184185
test_func(ElementJacobian2d ./ElementJacobian2d)
185186
test_func(CountIntegrator ./CountIntegrator)
186187
test_func(OmegahElementTests ./OmegahElementTests)
188+
test_func(PointMapping ./PointMapping)
187189
if(MeshFields_USE_EXCEPTIONS)
188190
# exception caught - no error
189191
test_func(ExceptionTest ./ExceptionTest)

‎src/MeshField.hpp‎

Lines changed: 8 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -81,13 +81,14 @@ struct LinearTriangleToVertexField {
8181
const auto vtxDim = 0;
8282
const auto ignored = -1;
8383
const auto localVtxIdx =
84-
Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored);
84+
(Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored) +
85+
2) %
86+
3;
8587
const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim);
8688
const MeshField::LO vtx = triVerts[(tri * triToVtxDegree) + localVtxIdx];
8789
return {0, triCompIdx, vtx, MeshField::Vertex}; // node, comp, ent, topo
8890
}
8991
};
90-
9192
struct QuadraticTriangleToField {
9293
Omega_h::LOs triVerts;
9394
Omega_h::LOs triEdges;
@@ -129,8 +130,10 @@ struct QuadraticTriangleToField {
129130
const auto triDim = 2;
130131
const auto vtxDim = 0;
131132
const auto ignored = -1;
132-
const auto localVtxIdx =
133-
Omega_h::simplex_down_template(triDim, vtxDim, dofHolderIdx, ignored);
133+
const auto localVtxIdx = (Omega_h::simplex_down_template(
134+
triDim, vtxDim, dofHolderIdx, ignored) +
135+
2) %
136+
3;
134137
const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim);
135138
osh_ent = triVerts[(tri * triToVtxDegree) + localVtxIdx];
136139
} else if (dofHolderTopo == MeshField::Edge) {
@@ -140,7 +143,7 @@ struct QuadraticTriangleToField {
140143
// passing dofHolderIdx as Omega_h_simplex.hpp does not provide
141144
// a function that maps a triangle and edge index to a 'canonical' edge
142145
// index. This may need to be revisited...
143-
osh_ent = triEdges[(tri * triToEdgeDegree) + dofHolderIdx];
146+
osh_ent = triEdges[(tri * triToEdgeDegree) + (dofHolderIdx + 2) % 3];
144147
} else {
145148
assert(false);
146149
}

‎test/testPointMapping.cpp‎

Lines changed: 122 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,122 @@
1+
#include <Kokkos_Core.hpp>
2+
#include <MeshField.hpp>
3+
#include <Omega_h_build.hpp>
4+
#include <Omega_h_mesh.hpp>
5+
template <size_t order> void runTest(Omega_h::Mesh mesh) {
6+
MeshField::OmegahMeshField<Kokkos::DefaultExecutionSpace, 2> mesh_field(mesh);
7+
auto shape_field =
8+
mesh_field.template CreateLagrangeField<double, order, 1>();
9+
10+
// initialize vertices
11+
auto dim = mesh.dim();
12+
auto coords = mesh.coords();
13+
auto f = KOKKOS_LAMBDA(double x, double y)->double { return 2 * x + 3 * y; };
14+
auto edge2vtx = mesh.ask_down(1, 0).ab2b;
15+
auto edgeMap = mesh.ask_down(dim, 1).ab2b;
16+
Kokkos::parallel_for(
17+
mesh.nverts(), KOKKOS_LAMBDA(int vtx) {
18+
auto x = coords[vtx * dim];
19+
auto y = coords[vtx * dim + 1];
20+
shape_field(vtx, 0, 0, MeshField::Mesh_Topology::Vertex) = f(x, y);
21+
});
22+
if (order == 2) {
23+
Kokkos::parallel_for(
24+
mesh.nedges(), KOKKOS_LAMBDA(int edge) {
25+
const auto left = edge2vtx[edge * 2];
26+
const auto right = edge2vtx[edge * 2 + 1];
27+
const auto x = (coords[left * dim] + coords[right * dim]) / 2.0;
28+
const auto y =
29+
(coords[left * dim + 1] + coords[right * dim + 1]) / 2.0;
30+
shape_field(edge, 0, 0, MeshField::Mesh_Topology::Edge) = f(x, y);
31+
});
32+
}
33+
const auto numNodesPerElem = order == 2 ? 6 : 3;
34+
Kokkos::View<double **> local_coords("", mesh.nelems() * numNodesPerElem, 3);
35+
Kokkos::parallel_for(
36+
"set", mesh.nelems() * numNodesPerElem, KOKKOS_LAMBDA(const int &i) {
37+
const auto val = i % numNodesPerElem;
38+
local_coords(i, 0) = (val == 0);
39+
local_coords(i, 1) = (val == 1);
40+
local_coords(i, 2) = (val == 2);
41+
if constexpr (order == 2) {
42+
if (val == 3) {
43+
local_coords(i, 0) = 1 / 2.0;
44+
local_coords(i, 1) = 1 / 2.0;
45+
} else if (val == 4) {
46+
local_coords(i, 1) = 1 / 2.0;
47+
local_coords(i, 2) = 1 / 2.0;
48+
} else if (val == 5) {
49+
local_coords(i, 0) = 1 / 2.0;
50+
local_coords(i, 2) = 1 / 2.0;
51+
}
52+
}
53+
});
54+
55+
auto eval_results = mesh_field.triangleLocalPointEval(
56+
local_coords, numNodesPerElem, shape_field);
57+
58+
int errors = 0;
59+
const auto triVerts = mesh.ask_elem_verts();
60+
const auto coordField = mesh_field.getCoordField();
61+
Kokkos::parallel_reduce(
62+
"test", mesh.nelems(),
63+
KOKKOS_LAMBDA(const int &tri, int &errors) {
64+
for (int node = 0; node < 3; ++node) {
65+
const auto triDim = 2;
66+
const auto vtxDim = 0;
67+
const auto ignored = -1;
68+
const auto localVtxIdx =
69+
Omega_h::simplex_down_template(triDim, vtxDim, node, ignored);
70+
const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim);
71+
int vtx = triVerts[(tri * triToVtxDegree) + localVtxIdx];
72+
auto x = coords[vtx * dim];
73+
auto y = coords[vtx * dim + 1];
74+
auto expected = f(x, y);
75+
auto actual = eval_results(tri * numNodesPerElem + node, 0);
76+
if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) {
77+
++errors;
78+
Kokkos::printf(
79+
"expected: %lf, actual: %lf, element: %d, node(vtx): %d\n",
80+
expected, actual, tri, node);
81+
}
82+
}
83+
for (int node = 3; node < numNodesPerElem; ++node) {
84+
const auto triDim = 2;
85+
const auto edgeDim = 1;
86+
const auto triToEdgeDegree = Omega_h::simplex_degree(triDim, edgeDim);
87+
const MeshField::LO triNode2DofHolder[3] = {0, 1, 2};
88+
int edge =
89+
edgeMap[(tri * triToEdgeDegree) + triNode2DofHolder[node - 3]];
90+
auto left = edge2vtx[edge * 2];
91+
auto right = edge2vtx[edge * 2 + 1];
92+
auto x = (coords[left * dim] + coords[right * dim]) / 2.0;
93+
auto y = (coords[left * dim + 1] + coords[right * dim + 1]) / 2.0;
94+
auto expected = f(x, y);
95+
auto actual = eval_results(tri * numNodesPerElem + node, 0);
96+
if (Kokkos::fabs(expected - actual) > MeshField::MachinePrecision) {
97+
++errors;
98+
Kokkos::printf(
99+
"expected: %lf, actual: %lf, element: %d, node(edge): %d\n",
100+
expected, actual, tri, node);
101+
}
102+
}
103+
},
104+
errors);
105+
if (errors > 0) {
106+
MeshField::fail("One or more mappings did not match\n");
107+
}
108+
}
109+
int main(int argc, char **argv) {
110+
// setup
111+
Kokkos::initialize(argc, argv);
112+
auto lib = Omega_h::Library(&argc, &argv);
113+
{
114+
auto world = lib.world();
115+
auto mesh =
116+
Omega_h::build_box(world, OMEGA_H_SIMPLEX, 1, 1, 1, 10, 10, 0, false);
117+
runTest<1>(mesh);
118+
runTest<2>(mesh);
119+
}
120+
Kokkos::finalize();
121+
return 0;
122+
}

0 commit comments

Comments
 (0)