From 326989561f7fccaed62ba2a900238b9bc9ce3c59 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Fri, 9 Oct 2026 14:46:47 -0400 Subject: [PATCH 01/12] Add the Wedge topology Add Wedge to Mesh_Topology and numWedge to MeshInfo for wedge (triangular prism) elements. Wedge is appended after Pyramid so no existing enum value changes; Prism and numPrism are unchanged. --- src/MeshField_Defines.hpp | 3 ++- src/MeshField_ShapeField.hpp | 1 + 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/src/MeshField_Defines.hpp b/src/MeshField_Defines.hpp index 8467a0b..1c79020 100644 --- a/src/MeshField_Defines.hpp +++ b/src/MeshField_Defines.hpp @@ -13,7 +13,8 @@ enum Mesh_Topology { Tetrahedron, Hexahedron, Prism, - Pyramid + Pyramid, + Wedge }; static bool Debug = false; const Real MachinePrecision = 1e-15; diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 5175b02..90ab0ea 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -25,6 +25,7 @@ struct MeshInfo { int numHex; // entDim = 3 int numPrism; // entDim = 3 int numPyramid; // entDim = 3 + int numWedge; // entDim = 3 int dim; }; From 5706a578797e2b3c6a76f16438253185762e532d Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Fri, 9 Oct 2026 14:46:47 -0400 Subject: [PATCH 02/12] Add LinearWedgeShape (P1 tensor-product wedge) LinearWedgeShape is the linear wedge (triangular prism): the tensor product of the linear triangle in (xi0, xi1) and the linear edge in zeta in [-1, 1]. Nodes 0-2 are the triangle vertices at zeta = -1 and nodes 3-5 those at zeta = +1, following Zienkiewicz, Taylor and Zhu, Section 6.4, Fig. 6.16a. test/testWedgeShape.cpp checks the Kronecker delta property at the nodes and, at interior points, partition of unity, reproduction of the bilinear functions spanned by the basis and their gradients, and the analytic gradients against finite differences. test/testWedgeElement.cpp builds a FieldElement from a coordinate field, LinearWedgeShape and a small wedge-to-vertex map, and checks MeshField::evaluate, getJacobians and getJacobianDeterminants on a right prism and a skewed wedge against the geometric map written out by hand. --- CMakeLists.txt | 4 + src/MeshField_Shape.hpp | 102 +++++++++++++++++++ test/testWedgeElement.cpp | 199 ++++++++++++++++++++++++++++++++++++++ test/testWedgeShape.cpp | 133 +++++++++++++++++++++++++ 4 files changed, 438 insertions(+) create mode 100644 test/testWedgeElement.cpp create mode 100644 test/testWedgeShape.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 45dc805..b60c19c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -213,6 +213,8 @@ meshfields_add_smoke_exe(OmegahTriTests test/testOmegahTri.cpp) meshfields_add_exe(ExceptionTest test/testExceptions.cpp) meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp) meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp) +meshfields_add_exe(WedgeShapeTest test/testWedgeShape.cpp) +meshfields_add_exe(WedgeElementTest test/testWedgeElement.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) if(MeshFields_USE_Cabana) @@ -235,6 +237,8 @@ test_func(CountIntegrator ./CountIntegrator) smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(OmegahTetTest ./OmegahTetTest) test_func(IntegratorMeasure ./IntegratorMeasure) +test_func(WedgeShapeTest ./WedgeShapeTest) +test_func(WedgeElementTest ./WedgeElementTest) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index 6fd4d47..ea977f1 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -557,5 +557,107 @@ struct QuadraticTetrahedronShape { } }; +/** + * @brief Linear (P1) shape functions for 3D wedge (triangular prism) elements + * + * Tensor product of the linear triangle in (xi0, xi1) and the linear edge in + * zeta = xi2. Defines 6 nodes, all located at the wedge vertices. + * Parametric coordinate range: xi0, xi1 \f$\in\f$ [0, 1] with + * xi0 + xi1 \f$\le\f$ 1, and zeta \f$\in\f$ [-1, 1] + * + * Node ordering: + * - Nodes 0-2: triangle vertices on the zeta = -1 face, at + * (xi0, xi1) = (0, 0), (1, 0), (0, 1) + * - Nodes 3-5: triangle vertices on the zeta = +1 face; node a+3 is joined to + * node a by the edge along zeta + * + * Shape functions, with L0 = 1-xi0-xi1, L1 = xi0, L2 = xi1 and a = 0, 1, 2: + * - \f$N_a = L_a (1 - \zeta) / 2\f$ + * - \f$N_{a+3} = L_a (1 + \zeta) / 2\f$ + * + * Shape functions and ordering from: + * Zienkiewicz, Taylor, and Zhu + * 'The Finite Element Method: Its Basis and Fundamentals', 2013, Section 6.4, + * Fig. 6.16a: \f$N_a = \frac{1}{2} L_a (1 + \zeta_a \zeta)\f$ + */ +struct LinearWedgeShape { + // clang-format off + static const size_t numNodes = 6; ///< Number of nodes (6 for linear wedge) + static const size_t meshEntDim = 3; ///< Mesh entity dimension (3D) + constexpr static Mesh_Topology DofHolders[1] = {Vertex}; ///< DOFs located at vertices + constexpr static size_t Order = 1; ///< Polynomial order (linear) + // clang-format on + + /** + * @brief Get parametric coordinates of element nodes + * @return Array of node coordinates [node0_xi0, node0_xi1, node0_zeta, ...] + */ + KOKKOS_INLINE_FUNCTION + Kokkos::Array getNodeParametricCoords() const { + // clang-format off + return {0,0,-1, //node 0 + 1,0,-1, //node 1 + 0,1,-1, //node 2 + 0,0, 1, //node 3 + 1,0, 1, //node 4 + 0,1, 1}; //node 5 + // clang-format on + } + + /** + * @brief Evaluate shape functions at parametric coordinate + * @param xi Parametric coordinates (xi0, xi1, zeta) where xi0, xi1 + * \f$\in\f$ [0,1], xi0 + xi1 \f$\le\f$ 1 and zeta \f$\in\f$ [-1,1] + * @return Array of shape function values [N0, ..., N5] + */ + KOKKOS_INLINE_FUNCTION + Kokkos::Array getValues(Vector3 const &xi) const { + assert(greaterThanOrEqual(xi[0], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[1], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[2], -1.0, ParametricCoordTol)); + assert(lessThanOrEqual(xi[2], 1.0, ParametricCoordTol)); + const Real L0 = 1 - xi[0] - xi[1]; + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + const Real down = (1 - xi[2]) / 2; // linear edge function of the -1 face + const Real up = (1 + xi[2]) / 2; // linear edge function of the +1 face + // clang-format off + return {L0*down, //node 0 + xi[0]*down, //node 1 + xi[1]*down, //node 2 + L0*up, //node 3 + xi[0]*up, //node 4 + xi[1]*up}; //node 5 + // clang-format on + } + + /** + * @brief Get gradients of shape functions in parametric coordinates + * @param xi Parametric coordinates (xi0, xi1, zeta), see getValues + * @return Array of gradients [dN0/dxi0, dN0/dxi1, dN0/dzeta, dN1/dxi0, ...] + */ + KOKKOS_INLINE_FUNCTION + Kokkos::Array + getLocalGradients(Vector3 const &xi) const { + assert(greaterThanOrEqual(xi[0], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[1], 0.0, ParametricCoordTol)); + assert(greaterThanOrEqual(xi[2], -1.0, ParametricCoordTol)); + assert(lessThanOrEqual(xi[2], 1.0, ParametricCoordTol)); + const Real L0 = 1 - xi[0] - xi[1]; + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + const Real L1 = xi[0]; + const Real L2 = xi[1]; + const Real down = (1 - xi[2]) / 2; + const Real up = (1 + xi[2]) / 2; + // clang-format off + return {-down, -down, -L0/2, //node 0 + down, 0, -L1/2, //node 1 + 0, down, -L2/2, //node 2 + -up, -up, L0/2, //node 3 + up, 0, L1/2, //node 4 + 0, up, L2/2}; //node 5 + // clang-format on + } +}; + } // namespace MeshField #endif diff --git a/test/testWedgeElement.cpp b/test/testWedgeElement.cpp new file mode 100644 index 0000000..245644c --- /dev/null +++ b/test/testWedgeElement.cpp @@ -0,0 +1,199 @@ +#include "KokkosController.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_Field.hpp" +#include "MeshField_For.hpp" +#include "MeshField_ShapeField.hpp" +#include +#include +#include +#include +#include + +// Checks MeshField::LinearWedgeShape inside a MeshField::FieldElement: field +// evaluation (MeshField::evaluate) and the Jacobian (MeshField::getJacobians, +// MeshField::getJacobianDeterminants) for two wedges, compared against the +// geometric map x(xi) differentiated by hand. Failures are counted explicitly +// so the test also fails in builds where assert() is disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; +using Vec3 = std::array; + +namespace { + +constexpr int numWedges = 2; +constexpr int numNodes = MeshField::LinearWedgeShape::numNodes; + +// Each wedge has its own six vertices: node n of wedge w is vertex 6*w + n. +// Nodes 0-2 are the zeta = -1 triangle and nodes 3-5 the zeta = +1 triangle, +// matching the LinearWedgeShape node ordering. +struct WedgeToVertexField { + static constexpr KOKKOS_FUNCTION Kokkos::Array + getTopology() { + return {MeshField::Wedge}; + } + + KOKKOS_FUNCTION MeshField::ElementToDofHolderMap + operator()(MeshField::LO wedgeNodeIdx, MeshField::LO wedgeCompIdx, + MeshField::LO wedge, MeshField::Mesh_Topology topo) const { + assert(topo == MeshField::Wedge); + const MeshField::LO vtx = numNodes * wedge + wedgeNodeIdx; + return {0, wedgeCompIdx, vtx, MeshField::Vertex}; + } +}; + +// clang-format off +// wedge 0: right prism, triangle (0,0),(2,0),(0,3) swept by height 4 +// x(xi) = (2 xi0, 3 xi1, 2 (1 + zeta)), J = diag(2, 3, 2), det J = 12 +// wedge 1: the top triangle is shifted, scaled and tilted, so J varies with +// xi0, xi1 and zeta +const Vec3 wedgeCoords[numWedges][numNodes] = { + {{0, 0, 0}, {2, 0, 0}, {0, 3, 0}, + {0, 0, 4}, {2, 0, 4}, {0, 3, 4}}, + {{0, 0, 0}, {1, 0, 0}, {0, 1, 0}, + {0.2, 0.1, 1}, {1.7, 0.3, 1.2}, {0.1, 1.4, 0.9}}}; +const char *wedgeNames[numWedges] = {"right prism", "skewed wedge"}; + +// the six nodes, the centroid, and interior points +const Vec3 localPts[] = { + {0, 0, -1}, {1, 0, -1}, {0, 1, -1}, {0, 0, 1}, {1, 0, 1}, {0, 1, 1}, + {1. / 3, 1. / 3, 0}, + {0.2, 0.3, -0.4}, {0.6, 0.1, 0.7}, {0.05, 0.9, 0.95}}; +// clang-format on +constexpr int numPts = sizeof(localPts) / sizeof(localPts[0]); +constexpr int centroidPt = 6; + +Vec3 add(Vec3 a, Vec3 b) { return {a[0] + b[0], a[1] + b[1], a[2] + b[2]}; } +Vec3 sub(Vec3 a, Vec3 b) { return {a[0] - b[0], a[1] - b[1], a[2] - b[2]}; } +Vec3 scale(Real s, Vec3 a) { return {s * a[0], s * a[1], s * a[2]}; } +Vec3 cross(Vec3 a, Vec3 b) { + return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], + a[0] * b[1] - a[1] * b[0]}; +} +Real dot(Vec3 a, Vec3 b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } + +// x(xi) = sum_a L_a(xi0, xi1) [ (1-zeta)/2 B_a + (1+zeta)/2 T_a ] with B_a the +// bottom (nodes 0-2) and T_a the top (nodes 3-5) vertex coordinates +Vec3 expectedX(const Vec3 (&nodes)[numNodes], Vec3 xi) { + const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; + Vec3 x = {0, 0, 0}; + for (int a = 0; a < 3; ++a) { + x = add(x, scale(L[a] * (1 - xi[2]) / 2, nodes[a])); + x = add(x, scale(L[a] * (1 + xi[2]) / 2, nodes[a + 3])); + } + return x; +} + +// rows of the MeshField Jacobian, J(i,j) = dx_j/dxi_i, from differentiating +// x(xi) directly: +// dx/dxi0 = (1-zeta)/2 (B1 - B0) + (1+zeta)/2 (T1 - T0) +// dx/dxi1 = (1-zeta)/2 (B2 - B0) + (1+zeta)/2 (T2 - T0) +// dx/dzeta = 1/2 sum_a L_a (T_a - B_a) +std::array expectedJ(const Vec3 (&nodes)[numNodes], Vec3 xi) { + const Real down = (1 - xi[2]) / 2; + const Real up = (1 + xi[2]) / 2; + const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; + const Vec3 dXi0 = add(scale(down, sub(nodes[1], nodes[0])), + scale(up, sub(nodes[4], nodes[3]))); + const Vec3 dXi1 = add(scale(down, sub(nodes[2], nodes[0])), + scale(up, sub(nodes[5], nodes[3]))); + Vec3 dZeta = {0, 0, 0}; + for (int a = 0; a < 3; ++a) + dZeta = add(dZeta, scale(L[a] / 2, sub(nodes[a + 3], nodes[a]))); + return {dXi0, dXi1, dZeta}; +} + +bool isClose(Real a, Real b) { + return std::fabs(a - b) <= 1e-12 * (1 + std::fabs(b)); +} + +template void setVtxCoords(CoordField field) { + const int numVtx = numWedges * numNodes; + Kokkos::View coords_h("coords_h", numVtx); + for (int w = 0; w < numWedges; ++w) + for (int n = 0; n < numNodes; ++n) + for (int d = 0; d < 3; ++d) + coords_h(numNodes * w + n, d) = wedgeCoords[w][n][d]; + auto coords = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), coords_h); + auto setCoordField = KOKKOS_LAMBDA(const int &vtx) { + for (int d = 0; d < 3; ++d) + field(vtx, 0, d, MeshField::Vertex) = coords(vtx, d); + }; + MeshField::parallel_for(ExecutionSpace(), {0}, {numVtx}, setCoordField, + "setCoordField"); +} + +} // namespace + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + int numFailed = 0; + { + MeshField::MeshInfo meshInfo{}; + meshInfo.numVtx = numWedges * numNodes; + meshInfo.numWedge = numWedges; + meshInfo.dim = 3; + auto coordField = MeshField::CreateCoordinateField< + ExecutionSpace, MeshField::KokkosController, 3>(meshInfo); + setVtxCoords(coordField.field); + + MeshField::FieldElement fe(numWedges, coordField.field, + MeshField::LinearWedgeShape(), + WedgeToVertexField()); + + // the same local points in every wedge, stored element-major + Kokkos::View lc("localCoords", numWedges * numPts, 3); + auto lc_h = Kokkos::create_mirror_view(lc); + for (int w = 0; w < numWedges; ++w) + for (int p = 0; p < numPts; ++p) + for (int d = 0; d < 3; ++d) + lc_h(w * numPts + p, d) = localPts[p][d]; + Kokkos::deep_copy(lc, lc_h); + + const auto x = MeshField::evaluate(fe, lc, numPts); + const auto J = MeshField::getJacobians(fe, lc, numPts); + const auto detJ = MeshField::getJacobianDeterminants(fe, J); + const auto x_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x); + const auto J_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), J); + const auto detJ_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ); + + for (int w = 0; w < numWedges; ++w) { + const auto &nodes = wedgeCoords[w]; + int evalFailed = 0, jacFailed = 0, detFailed = 0; + for (int p = 0; p < numPts; ++p) { + const int row = w * numPts + p; + // evaluate: x at a node is that node, x at the centroid is the mean + // of the six vertices, and x matches the map x(xi) everywhere + Vec3 xExpected = expectedX(nodes, localPts[p]); + if (p < numNodes) + xExpected = nodes[p]; + if (p == centroidPt) { + xExpected = {0, 0, 0}; + for (int n = 0; n < numNodes; ++n) + xExpected = add(xExpected, scale(1.0 / numNodes, nodes[n])); + } + for (int d = 0; d < 3; ++d) + evalFailed += !isClose(x_h(row, d), xExpected[d]); + // Jacobian rows and determinant + const auto jExpected = expectedJ(nodes, localPts[p]); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); + Real detExpected = dot(jExpected[0], cross(jExpected[1], jExpected[2])); + if (w == 0) + detExpected = 2.0 * 3.0 * 4.0 / 2.0; // right prism: a*b*h/2 = 12 + detFailed += !isClose(detJ_h(row), detExpected); + } + std::cout << "LinearWedgeShape FieldElement, " << wedgeNames[w] + << ": evaluate failures " << evalFailed + << ", Jacobian failures " << jacFailed << ", det J failures " + << detFailed << "\n"; + numFailed += evalFailed + jacFailed + detFailed; + } + } + Kokkos::finalize(); + return numFailed ? 1 : 0; +} diff --git a/test/testWedgeShape.cpp b/test/testWedgeShape.cpp new file mode 100644 index 0000000..c5804eb --- /dev/null +++ b/test/testWedgeShape.cpp @@ -0,0 +1,133 @@ +#include // before MeshField_Shape.hpp, which uses Kokkos +#include +#include +#include + +// Checks of MeshField::LinearWedgeShape. The checks run inside Kokkos kernels +// so the shape functions are exercised in the default execution space (e.g., +// on a GPU), and failures are counted explicitly so the test also fails in +// builds where assert() is disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; +using Shape = MeshField::LinearWedgeShape; + +namespace { + +KOKKOS_INLINE_FUNCTION bool isClose(Real a, Real b, Real tol) { + return Kokkos::fabs(a - b) <= tol; +} + +// A function in the span of the wedge basis, +// {1, xi0, xi1, zeta, xi0*zeta, xi1*zeta}, and its gradient +KOKKOS_INLINE_FUNCTION Real f(Real x0, Real x1, Real z) { + return 0.5 + 2 * x0 - 3 * x1 + 0.25 * z + 4 * x0 * z - 1.5 * x1 * z; +} +KOKKOS_INLINE_FUNCTION MeshField::Vector3 gradF(Real x0, Real x1, Real z) { + return {2 + 4 * z, -3 - 1.5 * z, 0.25 + 4 * x0 - 1.5 * x1}; +} + +// Each shape function is one at its own node and zero at the others. +int checkKroneckerDelta() { + int numFailed = 0; + Kokkos::parallel_reduce( + "wedgeKronecker", Kokkos::RangePolicy(0, Shape::numNodes), + KOKKOS_LAMBDA(const int node, int &failed) { + const Shape shp; + const auto coords = shp.getNodeParametricCoords(); + const MeshField::Vector3 xi{coords[3 * node], coords[3 * node + 1], + coords[3 * node + 2]}; + const auto N = shp.getValues(xi); + for (size_t m = 0; m < Shape::numNodes; ++m) { + const Real expected = (static_cast(node) == m) ? 1.0 : 0.0; + failed += !isClose(N[m], expected, MeshField::MachinePrecision); + } + }, + numFailed); + return numFailed; +} + +// At interior points: the values sum to one, the gradients sum to zero, the +// basis reproduces f and its gradient exactly, and the analytic gradients +// match central finite differences of the values. +int checkInteriorPoints() { + // clang-format off + const Real pts[][3] = {{1. / 3, 1. / 3, 0.0}, + {0.2, 0.3, -0.4}, + {0.6, 0.1, 0.7}, + {0.05, 0.9, 0.95}, + {0.45, 0.45, -0.95}}; + // clang-format on + const int numPts = sizeof(pts) / sizeof(pts[0]); + Kokkos::View pts_d("wedgePts", numPts); + auto pts_h = Kokkos::create_mirror_view(pts_d); + for (int i = 0; i < numPts; ++i) + for (int d = 0; d < 3; ++d) + pts_h(i, d) = pts[i][d]; + Kokkos::deep_copy(pts_d, pts_h); + + int numFailed = 0; + Kokkos::parallel_reduce( + "wedgeInterior", Kokkos::RangePolicy(0, numPts), + KOKKOS_LAMBDA(const int i, int &failed) { + const Shape shp; + const MeshField::Vector3 xi{pts_d(i, 0), pts_d(i, 1), pts_d(i, 2)}; + const auto N = shp.getValues(xi); + const auto G = shp.getLocalGradients(xi); + const auto nodes = shp.getNodeParametricCoords(); + + Real sumN = 0; + Real sumG[3] = {0, 0, 0}; + Real fInterp = 0; + Real gInterp[3] = {0, 0, 0}; + for (size_t n = 0; n < Shape::numNodes; ++n) { + const Real fn = f(nodes[3 * n], nodes[3 * n + 1], nodes[3 * n + 2]); + sumN += N[n]; + fInterp += N[n] * fn; + for (int d = 0; d < 3; ++d) { + sumG[d] += G[3 * n + d]; + gInterp[d] += G[3 * n + d] * fn; + } + } + const Real tol = 10 * MeshField::MachinePrecision; + failed += !isClose(sumN, 1.0, tol); + failed += !isClose(fInterp, f(xi[0], xi[1], xi[2]), tol); + const auto gExact = gradF(xi[0], xi[1], xi[2]); + for (int d = 0; d < 3; ++d) { + failed += !isClose(sumG[d], 0.0, tol); + failed += !isClose(gInterp[d], gExact[d], tol); + } + + const Real h = 1e-6; + for (int d = 0; d < 3; ++d) { + auto xp = xi; + auto xm = xi; + xp[d] += h; + xm[d] -= h; + const auto Np = shp.getValues(xp); + const auto Nm = shp.getValues(xm); + for (size_t n = 0; n < Shape::numNodes; ++n) + failed += !isClose((Np[n] - Nm[n]) / (2 * h), G[3 * n + d], 1e-8); + } + }, + numFailed); + return numFailed; +} + +} // namespace + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + int numFailed = 0; + { + const int kronecker = checkKroneckerDelta(); + std::cout << "LinearWedgeShape Kronecker delta failures: " << kronecker + << "\n"; + const int interior = checkInteriorPoints(); + std::cout << "LinearWedgeShape interior point failures: " << interior + << "\n"; + numFailed = kronecker + interior; + } + Kokkos::finalize(); + return numFailed ? 1 : 0; +} From 75b0f600b0f6519239e51f225014611d6f43c1e4 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 13 Oct 2026 11:51:54 -0400 Subject: [PATCH 03/12] improve LinearWedgeShape documentation --- src/MeshField_Shape.hpp | 28 +++++++++++++++++----------- 1 file changed, 17 insertions(+), 11 deletions(-) diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index ea977f1..8489814 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -558,27 +558,31 @@ struct QuadraticTetrahedronShape { }; /** - * @brief Linear (P1) shape functions for 3D wedge (triangular prism) elements + * @brief Linear shape functions for 3D wedge (triangular prism) elements * - * Tensor product of the linear triangle in (xi0, xi1) and the linear edge in - * zeta = xi2. Defines 6 nodes, all located at the wedge vertices. + * Defines linear basis functions for a wedge element with 6 nodes: the tensor product + * of the linear (P1) triangle in (xi0, xi1) and the linear edge in + * zeta = xi2. * Parametric coordinate range: xi0, xi1 \f$\in\f$ [0, 1] with * xi0 + xi1 \f$\le\f$ 1, and zeta \f$\in\f$ [-1, 1] * - * Node ordering: - * - Nodes 0-2: triangle vertices on the zeta = -1 face, at - * (xi0, xi1) = (0, 0), (1, 0), (0, 1) - * - Nodes 3-5: triangle vertices on the zeta = +1 face; node a+3 is joined to - * node a by the edge along zeta + * Node ordering (node a+3 is joined to node a by an edge along \zeta): + * - Node 0 at (0, 0 , -1) + * - Node 1 at (1, 0 , -1) + * - Node 2 at (0, 1 , -1) + * - Node 3 at (0, 0 , 1) + * - Node 4 at (1, 0 , 1) + * - Node 5 at (0, 1 , 1) * - * Shape functions, with L0 = 1-xi0-xi1, L1 = xi0, L2 = xi1 and a = 0, 1, 2: + * Shape functions use barycentric coordinates L0 = 1-xi0-xi1, L1 = xi0, L2 = xi1 of the triangle and, for a = 0, 1, 2: * - \f$N_a = L_a (1 - \zeta) / 2\f$ * - \f$N_{a+3} = L_a (1 + \zeta) / 2\f$ * * Shape functions and ordering from: * Zienkiewicz, Taylor, and Zhu * 'The Finite Element Method: Its Basis and Fundamentals', 2013, Section 6.4, - * Fig. 6.16a: \f$N_a = \frac{1}{2} L_a (1 + \zeta_a \zeta)\f$ + * Fig. 6.16a: \f$N_a = \frac{1}{2} L_a (1 + \zeta_a \zeta)\f$, with + * \f$\zeta_a = -1\f$ for nodes 0-2 and \f$+1\f$ for nodes 3-5 */ struct LinearWedgeShape { // clang-format off @@ -608,7 +612,9 @@ struct LinearWedgeShape { * @brief Evaluate shape functions at parametric coordinate * @param xi Parametric coordinates (xi0, xi1, zeta) where xi0, xi1 * \f$\in\f$ [0,1], xi0 + xi1 \f$\le\f$ 1 and zeta \f$\in\f$ [-1,1] - * @return Array of shape function values [N0, ..., N5] + * @return Array of shape function values [N0, ..., N5] where + * N_a = L_a ( 1 - \zeta)/2 and N_{a+3} = L_a (1+\zeta)/2 for a = 0,1,2 + * with L0 = 1 - x0 - xi1, L1 = xi0 and L2 = xi1 */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector3 const &xi) const { From 554a440f1af59260fdc5f7b9a58cd1dd264c425a Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 13 Oct 2026 12:03:59 -0400 Subject: [PATCH 04/12] refactor --- test/testWedgeElement.cpp | 5 +---- test/testWedgeShape.cpp | 2 +- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/test/testWedgeElement.cpp b/test/testWedgeElement.cpp index 245644c..9768a8a 100644 --- a/test/testWedgeElement.cpp +++ b/test/testWedgeElement.cpp @@ -129,10 +129,7 @@ int main(int argc, char **argv) { Kokkos::initialize(argc, argv); int numFailed = 0; { - MeshField::MeshInfo meshInfo{}; - meshInfo.numVtx = numWedges * numNodes; - meshInfo.numWedge = numWedges; - meshInfo.dim = 3; + const MeshField::MeshInfo meshInfo{.numVtx = numWedges * numNodes, .numWedge = numWedges, .dim = 3}; auto coordField = MeshField::CreateCoordinateField< ExecutionSpace, MeshField::KokkosController, 3>(meshInfo); setVtxCoords(coordField.field); diff --git a/test/testWedgeShape.cpp b/test/testWedgeShape.cpp index c5804eb..c946bb5 100644 --- a/test/testWedgeShape.cpp +++ b/test/testWedgeShape.cpp @@ -1,4 +1,4 @@ -#include // before MeshField_Shape.hpp, which uses Kokkos +#include #include #include #include From 6f3cf6f55b977032f2132b6c3f992a5753ed0746 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 13 Oct 2026 12:08:50 -0400 Subject: [PATCH 05/12] remove namespace --- test/testWedgeShape.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/testWedgeShape.cpp b/test/testWedgeShape.cpp index c946bb5..06ed98d 100644 --- a/test/testWedgeShape.cpp +++ b/test/testWedgeShape.cpp @@ -12,7 +12,6 @@ using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; using Shape = MeshField::LinearWedgeShape; -namespace { KOKKOS_INLINE_FUNCTION bool isClose(Real a, Real b, Real tol) { return Kokkos::fabs(a - b) <= tol; @@ -114,7 +113,6 @@ int checkInteriorPoints() { return numFailed; } -} // namespace int main(int argc, char **argv) { Kokkos::initialize(argc, argv); From 8d3550578057c46dd818b61c6475e3362ae1dedb Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 13 Oct 2026 16:12:11 -0400 Subject: [PATCH 06/12] replace custom vector functions with omega_h_vector functions --- test/testWedgeElement.cpp | 52 ++++++++++++++++----------------------- test/testWedgeShape.cpp | 4 +-- 2 files changed, 23 insertions(+), 33 deletions(-) diff --git a/test/testWedgeElement.cpp b/test/testWedgeElement.cpp index 9768a8a..9f15c43 100644 --- a/test/testWedgeElement.cpp +++ b/test/testWedgeElement.cpp @@ -3,11 +3,10 @@ #include "MeshField_Field.hpp" #include "MeshField_For.hpp" #include "MeshField_ShapeField.hpp" +#include #include #include -#include #include -#include // Checks MeshField::LinearWedgeShape inside a MeshField::FieldElement: field // evaluation (MeshField::evaluate) and the Jacobian (MeshField::getJacobians, @@ -17,13 +16,17 @@ using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; -using Vec3 = std::array; +using Vec3 = Omega_h::Vector<3>; -namespace { constexpr int numWedges = 2; constexpr int numNodes = MeshField::LinearWedgeShape::numNodes; +bool isClose(Real a, Real b) { + return Omega_h::are_close(a, b, 1e-12, 1e-12); +} + + // Each wedge has its own six vertices: node n of wedge w is vertex 6*w + n. // Nodes 0-2 are the zeta = -1 triangle and nodes 3-5 the zeta = +1 triangle, // matching the LinearWedgeShape node ordering. @@ -63,23 +66,14 @@ const Vec3 localPts[] = { constexpr int numPts = sizeof(localPts) / sizeof(localPts[0]); constexpr int centroidPt = 6; -Vec3 add(Vec3 a, Vec3 b) { return {a[0] + b[0], a[1] + b[1], a[2] + b[2]}; } -Vec3 sub(Vec3 a, Vec3 b) { return {a[0] - b[0], a[1] - b[1], a[2] - b[2]}; } -Vec3 scale(Real s, Vec3 a) { return {s * a[0], s * a[1], s * a[2]}; } -Vec3 cross(Vec3 a, Vec3 b) { - return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], - a[0] * b[1] - a[1] * b[0]}; -} -Real dot(Vec3 a, Vec3 b) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; } - // x(xi) = sum_a L_a(xi0, xi1) [ (1-zeta)/2 B_a + (1+zeta)/2 T_a ] with B_a the // bottom (nodes 0-2) and T_a the top (nodes 3-5) vertex coordinates Vec3 expectedX(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; - Vec3 x = {0, 0, 0}; + Vec3 x = Omega_h::zero_vector<3>(); for (int a = 0; a < 3; ++a) { - x = add(x, scale(L[a] * (1 - xi[2]) / 2, nodes[a])); - x = add(x, scale(L[a] * (1 + xi[2]) / 2, nodes[a + 3])); + x += L[a] * (1 - xi[2]) / 2 * nodes[a]; + x += L[a] * (1 + xi[2]) / 2 * nodes[a + 3]; } return x; } @@ -93,21 +87,18 @@ std::array expectedJ(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Real down = (1 - xi[2]) / 2; const Real up = (1 + xi[2]) / 2; const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; - const Vec3 dXi0 = add(scale(down, sub(nodes[1], nodes[0])), - scale(up, sub(nodes[4], nodes[3]))); - const Vec3 dXi1 = add(scale(down, sub(nodes[2], nodes[0])), - scale(up, sub(nodes[5], nodes[3]))); - Vec3 dZeta = {0, 0, 0}; - for (int a = 0; a < 3; ++a) - dZeta = add(dZeta, scale(L[a] / 2, sub(nodes[a + 3], nodes[a]))); + const Vec3 dXi0 = down * (nodes[1] - nodes[0]) + up * (nodes[4] - nodes[3]); + const Vec3 dXi1 = down * (nodes[2] - nodes[0]) + up * (nodes[5] - nodes[3]); + Vec3 dZeta = Omega_h::zero_vector<3>(); + for (int a = 0; a < 3; ++a){ + dZeta += L[a] / 2 * (nodes[a + 3] - nodes[a]); + } return {dXi0, dXi1, dZeta}; } -bool isClose(Real a, Real b) { - return std::fabs(a - b) <= 1e-12 * (1 + std::fabs(b)); -} -template void setVtxCoords(CoordField field) { +template +void setVtxCoords(CoordField field) { const int numVtx = numWedges * numNodes; Kokkos::View coords_h("coords_h", numVtx); for (int w = 0; w < numWedges; ++w) @@ -123,7 +114,6 @@ template void setVtxCoords(CoordField field) { "setCoordField"); } -} // namespace int main(int argc, char **argv) { Kokkos::initialize(argc, argv); @@ -168,9 +158,9 @@ int main(int argc, char **argv) { if (p < numNodes) xExpected = nodes[p]; if (p == centroidPt) { - xExpected = {0, 0, 0}; + xExpected = Omega_h::zero_vector<3>(); for (int n = 0; n < numNodes; ++n) - xExpected = add(xExpected, scale(1.0 / numNodes, nodes[n])); + xExpected += nodes[n]/numNodes ; } for (int d = 0; d < 3; ++d) evalFailed += !isClose(x_h(row, d), xExpected[d]); @@ -179,7 +169,7 @@ int main(int argc, char **argv) { for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); - Real detExpected = dot(jExpected[0], cross(jExpected[1], jExpected[2])); + Real detExpected = Omega_h::inner_product(jExpected[0] * Omega_h:: cross(jExpected[1], jExpected[2])); if (w == 0) detExpected = 2.0 * 3.0 * 4.0 / 2.0; // right prism: a*b*h/2 = 12 detFailed += !isClose(detJ_h(row), detExpected); diff --git a/test/testWedgeShape.cpp b/test/testWedgeShape.cpp index 06ed98d..4a12b6b 100644 --- a/test/testWedgeShape.cpp +++ b/test/testWedgeShape.cpp @@ -4,8 +4,8 @@ #include // Checks of MeshField::LinearWedgeShape. The checks run inside Kokkos kernels -// so the shape functions are exercised in the default execution space (e.g., -// on a GPU), and failures are counted explicitly so the test also fails in +// so the shape functions are exercised in the default execution space +// and failures are counted explicitly so the test also fails in // builds where assert() is disabled. using ExecutionSpace = Kokkos::DefaultExecutionSpace; From fffde66f92a7941eac42986bcd8b738d1afb2fed Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 14 Oct 2026 10:47:15 -0400 Subject: [PATCH 07/12] use omega_h_vector functions --- test/testWedgeElement.cpp | 2 +- test/testWedgeShape.cpp | 2 ++ 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/test/testWedgeElement.cpp b/test/testWedgeElement.cpp index 9f15c43..0731769 100644 --- a/test/testWedgeElement.cpp +++ b/test/testWedgeElement.cpp @@ -68,7 +68,7 @@ constexpr int centroidPt = 6; // x(xi) = sum_a L_a(xi0, xi1) [ (1-zeta)/2 B_a + (1+zeta)/2 T_a ] with B_a the // bottom (nodes 0-2) and T_a the top (nodes 3-5) vertex coordinates -Vec3 expectedX(const Vec3 (&nodes)[numNodes], Vec3 xi) { +Vec3 expectedX(const std::array& nodes, Vec3 xi) { const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; Vec3 x = Omega_h::zero_vector<3>(); for (int a = 0; a < 3; ++a) { diff --git a/test/testWedgeShape.cpp b/test/testWedgeShape.cpp index 4a12b6b..b5a8aa6 100644 --- a/test/testWedgeShape.cpp +++ b/test/testWedgeShape.cpp @@ -97,6 +97,8 @@ int checkInteriorPoints() { failed += !isClose(gInterp[d], gExact[d], tol); } + // gradient check with finite difference evaluation + const Real h = 1e-6; for (int d = 0; d < 3; ++d) { auto xp = xi; From 296d2beec0f7bd0e4a7ffe701dfc916d6bde18fb Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 14 Oct 2026 10:54:41 -0400 Subject: [PATCH 08/12] change to std::array --- test/testWedgeElement.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/testWedgeElement.cpp b/test/testWedgeElement.cpp index 0731769..ac89d46 100644 --- a/test/testWedgeElement.cpp +++ b/test/testWedgeElement.cpp @@ -76,6 +76,7 @@ Vec3 expectedX(const std::array& nodes, Vec3 xi) { x += L[a] * (1 + xi[2]) / 2 * nodes[a + 3]; } return x; + } // rows of the MeshField Jacobian, J(i,j) = dx_j/dxi_i, from differentiating @@ -83,7 +84,7 @@ Vec3 expectedX(const std::array& nodes, Vec3 xi) { // dx/dxi0 = (1-zeta)/2 (B1 - B0) + (1+zeta)/2 (T1 - T0) // dx/dxi1 = (1-zeta)/2 (B2 - B0) + (1+zeta)/2 (T2 - T0) // dx/dzeta = 1/2 sum_a L_a (T_a - B_a) -std::array expectedJ(const Vec3 (&nodes)[numNodes], Vec3 xi) { +std::array expectedJ(const std::array& nodes, Vec3 xi) { const Real down = (1 - xi[2]) / 2; const Real up = (1 + xi[2]) / 2; const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; @@ -176,6 +177,7 @@ int main(int argc, char **argv) { } std::cout << "LinearWedgeShape FieldElement, " << wedgeNames[w] << ": evaluate failures " << evalFailed + << ", Jacobian failures " << jacFailed << ", det J failures " << detFailed << "\n"; numFailed += evalFailed + jacFailed + detFailed; From b58f94f9d5bd3589f0c8fcbe555932ab062ebd77 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 14 Oct 2026 11:51:22 -0400 Subject: [PATCH 09/12] change to C-type array --- test/testWedgeElement.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/testWedgeElement.cpp b/test/testWedgeElement.cpp index ac89d46..405f3a2 100644 --- a/test/testWedgeElement.cpp +++ b/test/testWedgeElement.cpp @@ -68,7 +68,7 @@ constexpr int centroidPt = 6; // x(xi) = sum_a L_a(xi0, xi1) [ (1-zeta)/2 B_a + (1+zeta)/2 T_a ] with B_a the // bottom (nodes 0-2) and T_a the top (nodes 3-5) vertex coordinates -Vec3 expectedX(const std::array& nodes, Vec3 xi) { +Vec3 expectedX(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; Vec3 x = Omega_h::zero_vector<3>(); for (int a = 0; a < 3; ++a) { @@ -84,7 +84,7 @@ Vec3 expectedX(const std::array& nodes, Vec3 xi) { // dx/dxi0 = (1-zeta)/2 (B1 - B0) + (1+zeta)/2 (T1 - T0) // dx/dxi1 = (1-zeta)/2 (B2 - B0) + (1+zeta)/2 (T2 - T0) // dx/dzeta = 1/2 sum_a L_a (T_a - B_a) -std::array expectedJ(const std::array& nodes, Vec3 xi) { +std::array expectedJ(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Real down = (1 - xi[2]) / 2; const Real up = (1 + xi[2]) / 2; const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; @@ -170,7 +170,7 @@ int main(int argc, char **argv) { for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); - Real detExpected = Omega_h::inner_product(jExpected[0] * Omega_h:: cross(jExpected[1], jExpected[2])); + Real detExpected = jExpected[0] * Omega_h:: cross(jExpected[1], jExpected[2]); if (w == 0) detExpected = 2.0 * 3.0 * 4.0 / 2.0; // right prism: a*b*h/2 = 12 detFailed += !isClose(detJ_h(row), detExpected); From 54d471a0c64464a09a1b73fe3f3220df16dfc185 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 15 Oct 2026 14:12:05 -0400 Subject: [PATCH 10/12] rename wedge tests to match class names and split checks into run functions --- ...Element.cpp => testLinearWedgeElement.cpp} | 163 +++++++++--------- ...edgeShape.cpp => testLinearWedgeShape.cpp} | 53 +++--- 2 files changed, 114 insertions(+), 102 deletions(-) rename test/{testWedgeElement.cpp => testLinearWedgeElement.cpp} (54%) rename test/{testWedgeShape.cpp => testLinearWedgeShape.cpp} (80%) diff --git a/test/testWedgeElement.cpp b/test/testLinearWedgeElement.cpp similarity index 54% rename from test/testWedgeElement.cpp rename to test/testLinearWedgeElement.cpp index 405f3a2..0d46d3f 100644 --- a/test/testWedgeElement.cpp +++ b/test/testLinearWedgeElement.cpp @@ -3,9 +3,10 @@ #include "MeshField_Field.hpp" #include "MeshField_For.hpp" #include "MeshField_ShapeField.hpp" -#include #include +#include #include +#include #include // Checks MeshField::LinearWedgeShape inside a MeshField::FieldElement: field @@ -18,15 +19,13 @@ using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; using Vec3 = Omega_h::Vector<3>; - constexpr int numWedges = 2; constexpr int numNodes = MeshField::LinearWedgeShape::numNodes; -bool isClose(Real a, Real b) { - return Omega_h::are_close(a, b, 1e-12, 1e-12); +bool isClose(Real a, Real b) { + return Omega_h::are_close(a, b, 1e-12, 1e-12); } - // Each wedge has its own six vertices: node n of wedge w is vertex 6*w + n. // Nodes 0-2 are the zeta = -1 triangle and nodes 3-5 the zeta = +1 triangle, // matching the LinearWedgeShape node ordering. @@ -63,7 +62,7 @@ const Vec3 localPts[] = { {1. / 3, 1. / 3, 0}, {0.2, 0.3, -0.4}, {0.6, 0.1, 0.7}, {0.05, 0.9, 0.95}}; // clang-format on -constexpr int numPts = sizeof(localPts) / sizeof(localPts[0]); +constexpr int numPts = static_cast(sizeof(localPts) / sizeof(localPts[0])); constexpr int centroidPt = 6; // x(xi) = sum_a L_a(xi0, xi1) [ (1-zeta)/2 B_a + (1+zeta)/2 T_a ] with B_a the @@ -72,11 +71,10 @@ Vec3 expectedX(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; Vec3 x = Omega_h::zero_vector<3>(); for (int a = 0; a < 3; ++a) { - x += L[a] * (1 - xi[2]) / 2 * nodes[a]; - x += L[a] * (1 + xi[2]) / 2 * nodes[a + 3]; + x += L[a] * (1 - xi[2]) / 2 * nodes[a]; + x += L[a] * (1 + xi[2]) / 2 * nodes[a + 3]; } return x; - } // rows of the MeshField Jacobian, J(i,j) = dx_j/dxi_i, from differentiating @@ -91,14 +89,13 @@ std::array expectedJ(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Vec3 dXi0 = down * (nodes[1] - nodes[0]) + up * (nodes[4] - nodes[3]); const Vec3 dXi1 = down * (nodes[2] - nodes[0]) + up * (nodes[5] - nodes[3]); Vec3 dZeta = Omega_h::zero_vector<3>(); - for (int a = 0; a < 3; ++a){ + for (int a = 0; a < 3; ++a) dZeta += L[a] / 2 * (nodes[a + 3] - nodes[a]); - } return {dXi0, dXi1, dZeta}; } - -template +// Writes the wedge vertex coordinates into the coordinate field. +template void setVtxCoords(CoordField field) { const int numVtx = numWedges * numNodes; Kokkos::View coords_h("coords_h", numVtx); @@ -115,74 +112,82 @@ void setVtxCoords(CoordField field) { "setCoordField"); } +// Compares the evaluated positions, Jacobians and determinants of wedge w with +// the hand-differentiated map at every local point. Returns the failure count. +template +int checkWedge(int w, const XView &x_h, const JView &J_h, + const DetView &detJ_h) { + const auto &nodes = wedgeCoords[w]; + int evalFailed = 0, jacFailed = 0, detFailed = 0; + for (int p = 0; p < numPts; ++p) { + const int row = w * numPts + p; + // at a node, x is that node; at the centroid, x is the mean of the six + // vertices; everywhere else x follows the map + Vec3 xExpected = expectedX(nodes, localPts[p]); + if (p < numNodes) + xExpected = nodes[p]; + if (p == centroidPt) { + xExpected = Omega_h::zero_vector<3>(); + for (int n = 0; n < numNodes; ++n) + xExpected += nodes[n] / Real(numNodes); + } + for (int d = 0; d < 3; ++d) + evalFailed += !isClose(x_h(row, d), xExpected[d]); + + const auto jExpected = expectedJ(nodes, localPts[p]); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); + const Real detExpected = + jExpected[0] * Omega_h::cross(jExpected[1], jExpected[2]); + detFailed += !isClose(detJ_h(row), detExpected); + } + std::cout << "LinearWedgeShape FieldElement, " << wedgeNames[w] + << ": evaluate failures " << evalFailed << ", Jacobian failures " + << jacFailed << ", det J failures " << detFailed << "\n"; + return evalFailed + jacFailed + detFailed; +} + +// Builds the coordinate field and the FieldElement, runs the evaluation, +// Jacobian and determinant, and checks every wedge. Must be called between +// Kokkos::initialize and Kokkos::finalize. +int runLinearWedgeElementChecks() { + const MeshField::MeshInfo meshInfo{ + .numVtx = numWedges * numNodes, .numWedge = numWedges, .dim = 3}; + auto coordField = MeshField::CreateCoordinateField< + ExecutionSpace, MeshField::KokkosController, 3>(meshInfo); + setVtxCoords(coordField.field); + + MeshField::FieldElement fe(numWedges, coordField.field, + MeshField::LinearWedgeShape(), + WedgeToVertexField()); + + // the same local points in every wedge, stored element-major + Kokkos::View lc("localCoords", numWedges * numPts, 3); + auto lc_h = Kokkos::create_mirror_view(lc); + for (int w = 0; w < numWedges; ++w) + for (int p = 0; p < numPts; ++p) + for (int d = 0; d < 3; ++d) + lc_h(w * numPts + p, d) = localPts[p][d]; + Kokkos::deep_copy(lc, lc_h); + + const auto x = MeshField::evaluate(fe, lc, numPts); + const auto J = MeshField::getJacobians(fe, lc, numPts); + const auto detJ = MeshField::getJacobianDeterminants(fe, J); + const auto x_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x); + const auto J_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), J); + const auto detJ_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ); + + int numFailed = 0; + for (int w = 0; w < numWedges; ++w) + numFailed += checkWedge(w, x_h, J_h, detJ_h); + return numFailed; +} int main(int argc, char **argv) { Kokkos::initialize(argc, argv); - int numFailed = 0; - { - const MeshField::MeshInfo meshInfo{.numVtx = numWedges * numNodes, .numWedge = numWedges, .dim = 3}; - auto coordField = MeshField::CreateCoordinateField< - ExecutionSpace, MeshField::KokkosController, 3>(meshInfo); - setVtxCoords(coordField.field); - - MeshField::FieldElement fe(numWedges, coordField.field, - MeshField::LinearWedgeShape(), - WedgeToVertexField()); - - // the same local points in every wedge, stored element-major - Kokkos::View lc("localCoords", numWedges * numPts, 3); - auto lc_h = Kokkos::create_mirror_view(lc); - for (int w = 0; w < numWedges; ++w) - for (int p = 0; p < numPts; ++p) - for (int d = 0; d < 3; ++d) - lc_h(w * numPts + p, d) = localPts[p][d]; - Kokkos::deep_copy(lc, lc_h); - - const auto x = MeshField::evaluate(fe, lc, numPts); - const auto J = MeshField::getJacobians(fe, lc, numPts); - const auto detJ = MeshField::getJacobianDeterminants(fe, J); - const auto x_h = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x); - const auto J_h = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), J); - const auto detJ_h = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ); - - for (int w = 0; w < numWedges; ++w) { - const auto &nodes = wedgeCoords[w]; - int evalFailed = 0, jacFailed = 0, detFailed = 0; - for (int p = 0; p < numPts; ++p) { - const int row = w * numPts + p; - // evaluate: x at a node is that node, x at the centroid is the mean - // of the six vertices, and x matches the map x(xi) everywhere - Vec3 xExpected = expectedX(nodes, localPts[p]); - if (p < numNodes) - xExpected = nodes[p]; - if (p == centroidPt) { - xExpected = Omega_h::zero_vector<3>(); - for (int n = 0; n < numNodes; ++n) - xExpected += nodes[n]/numNodes ; - } - for (int d = 0; d < 3; ++d) - evalFailed += !isClose(x_h(row, d), xExpected[d]); - // Jacobian rows and determinant - const auto jExpected = expectedJ(nodes, localPts[p]); - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); - Real detExpected = jExpected[0] * Omega_h:: cross(jExpected[1], jExpected[2]); - if (w == 0) - detExpected = 2.0 * 3.0 * 4.0 / 2.0; // right prism: a*b*h/2 = 12 - detFailed += !isClose(detJ_h(row), detExpected); - } - std::cout << "LinearWedgeShape FieldElement, " << wedgeNames[w] - << ": evaluate failures " << evalFailed - - << ", Jacobian failures " << jacFailed << ", det J failures " - << detFailed << "\n"; - numFailed += evalFailed + jacFailed + detFailed; - } - } + const int numFailed = runLinearWedgeElementChecks(); Kokkos::finalize(); return numFailed ? 1 : 0; -} +} \ No newline at end of file diff --git a/test/testWedgeShape.cpp b/test/testLinearWedgeShape.cpp similarity index 80% rename from test/testWedgeShape.cpp rename to test/testLinearWedgeShape.cpp index b5a8aa6..267c966 100644 --- a/test/testWedgeShape.cpp +++ b/test/testLinearWedgeShape.cpp @@ -1,18 +1,17 @@ -#include +#include #include #include #include // Checks of MeshField::LinearWedgeShape. The checks run inside Kokkos kernels -// so the shape functions are exercised in the default execution space -// and failures are counted explicitly so the test also fails in -// builds where assert() is disabled. +// so the shape functions are exercised in the default execution space, and +// failures are counted explicitly so the test also fails in builds where +// assert() is disabled. using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; using Shape = MeshField::LinearWedgeShape; - KOKKOS_INLINE_FUNCTION bool isClose(Real a, Real b, Real tol) { return Kokkos::fabs(a - b) <= tol; } @@ -46,10 +45,8 @@ int checkKroneckerDelta() { return numFailed; } -// At interior points: the values sum to one, the gradients sum to zero, the -// basis reproduces f and its gradient exactly, and the analytic gradients -// match central finite differences of the values. -int checkInteriorPoints() { +// Interior points, spanning the zeta range. +Kokkos::View makeInteriorPoints() { // clang-format off const Real pts[][3] = {{1. / 3, 1. / 3, 0.0}, {0.2, 0.3, -0.4}, @@ -57,13 +54,22 @@ int checkInteriorPoints() { {0.05, 0.9, 0.95}, {0.45, 0.45, -0.95}}; // clang-format on - const int numPts = sizeof(pts) / sizeof(pts[0]); + const int numPts = static_cast(sizeof(pts) / sizeof(pts[0])); Kokkos::View pts_d("wedgePts", numPts); auto pts_h = Kokkos::create_mirror_view(pts_d); for (int i = 0; i < numPts; ++i) for (int d = 0; d < 3; ++d) pts_h(i, d) = pts[i][d]; Kokkos::deep_copy(pts_d, pts_h); + return pts_d; +} + +// At interior points: the values sum to one, the gradients sum to zero, the +// basis reproduces f and its gradient exactly, and the analytic gradients +// match central finite differences of the values. +int checkInteriorPoints() { + const auto pts_d = makeInteriorPoints(); + const int numPts = static_cast(pts_d.extent(0)); int numFailed = 0; Kokkos::parallel_reduce( @@ -97,8 +103,7 @@ int checkInteriorPoints() { failed += !isClose(gInterp[d], gExact[d], tol); } - // gradient check with finite difference evaluation - + // gradients against central finite differences of the values const Real h = 1e-6; for (int d = 0; d < 3; ++d) { auto xp = xi; @@ -115,19 +120,21 @@ int checkInteriorPoints() { return numFailed; } +// Runs both checks and prints their failure counts. Must be called between +// Kokkos::initialize and Kokkos::finalize. +int runLinearWedgeShapeChecks() { + const int kronecker = checkKroneckerDelta(); + std::cout << "LinearWedgeShape Kronecker delta failures: " << kronecker + << "\n"; + const int interior = checkInteriorPoints(); + std::cout << "LinearWedgeShape interior point failures: " << interior + << "\n"; + return kronecker + interior; +} int main(int argc, char **argv) { Kokkos::initialize(argc, argv); - int numFailed = 0; - { - const int kronecker = checkKroneckerDelta(); - std::cout << "LinearWedgeShape Kronecker delta failures: " << kronecker - << "\n"; - const int interior = checkInteriorPoints(); - std::cout << "LinearWedgeShape interior point failures: " << interior - << "\n"; - numFailed = kronecker + interior; - } + const int numFailed = runLinearWedgeShapeChecks(); Kokkos::finalize(); return numFailed ? 1 : 0; -} +} \ No newline at end of file From a198061e36a77877dd6734d2fc946a1d40a96ab0 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 15 Oct 2026 14:18:15 -0400 Subject: [PATCH 11/12] update the name change in CMakeLists --- CMakeLists.txt | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b60c19c..3376e28 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -213,8 +213,8 @@ meshfields_add_smoke_exe(OmegahTriTests test/testOmegahTri.cpp) meshfields_add_exe(ExceptionTest test/testExceptions.cpp) meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp) meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp) -meshfields_add_exe(WedgeShapeTest test/testWedgeShape.cpp) -meshfields_add_exe(WedgeElementTest test/testWedgeElement.cpp) +meshfields_add_exe(LinearWedgeShapeTest test/testLinearWedgeShape.cpp) +meshfields_add_exe(LinearWedgeElementTest test/testLinearWedgeElement.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) if(MeshFields_USE_Cabana) @@ -237,8 +237,8 @@ test_func(CountIntegrator ./CountIntegrator) smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(OmegahTetTest ./OmegahTetTest) test_func(IntegratorMeasure ./IntegratorMeasure) -test_func(WedgeShapeTest ./WedgeShapeTest) -test_func(WedgeElementTest ./WedgeElementTest) +test_func(LinearWedgeShapeTest ./LinearWedgeShapeTest) +test_func(LinearWedgeElementTest ./LinearWedgeElementTest) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) From 86715bf888b4964832e5ce4e7067eef2be15ae33 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 28 Oct 2026 11:05:33 -0400 Subject: [PATCH 12/12] Use the controller-based MeshField::parallel_for in LinearWedgeElementTest --- test/testLinearWedgeElement.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/testLinearWedgeElement.cpp b/test/testLinearWedgeElement.cpp index 0d46d3f..057a51a 100644 --- a/test/testLinearWedgeElement.cpp +++ b/test/testLinearWedgeElement.cpp @@ -108,8 +108,8 @@ void setVtxCoords(CoordField field) { for (int d = 0; d < 3; ++d) field(vtx, 0, d, MeshField::Vertex) = coords(vtx, d); }; - MeshField::parallel_for(ExecutionSpace(), {0}, {numVtx}, setCoordField, - "setCoordField"); + MeshField::parallel_for( + {0}, {numVtx}, setCoordField, "setCoordField"); } // Compares the evaluated positions, Jacobians and determinants of wedge w with