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/23] 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/23] 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/23] 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/23] 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/23] 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/23] 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/23] 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/23] 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/23] 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 f4db17e915b74e6a79be156c85853b038121a8e9 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 14 Oct 2026 15:15:22 -0400 Subject: [PATCH 10/23] Add LinearCubicWedgeShape for cubic wedge elements Add LinearCubicWedgeShape, a 12-node wedge with a linear triangle basis in (xi0, xi1) and a cubic Lagrange basis in zeta. The shape is the tensor product of the two. Nodes 0-5 match LinearWedgeShape: triangle vertices at zeta = -1 and +1. Nodes 6-11 are two interior nodes on each vertical edge, at zeta = -1/3 and +1/3, so each edge carries two DOFs. The struct declares DofHolders {Vertex, Edge}, NumDofHolders {6, 3}, DofsPerHolder {1, 2}, and Order = 3. The cubic functions are the Lagrange interpolants on the stations -1, -1/3, 1/3, 1. Gradients use their analytic derivatives. --- src/MeshField_Shape.hpp | 163 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 163 insertions(+) diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index 8489814..419c70e 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -665,5 +665,168 @@ struct LinearWedgeShape { } }; + +/** + * @brief Linear in the triangle and cubic in the zeta direction shape functions for 3D tensor product wedge elements + * + * Defines shape functions for a wedge element with 12 nodes: the tensor product + * of the linear triangle basis in (xi0, xi1) and the cubic lagrange basis in + * zeta = xi2. + * + * Parametric coordinate range: xi0, xi1 \f$\in\f$ [0, 1] with + * xi0, xi1 \f$\ge\f$ 0, xi0 + xi1 \f$\le\f$ 1, and zeta \f$\in\f$ [-1, 1] + * + * Node ordering (nodes 0-5 as in LinearWedgeShape): + * - Node 0-2: triangle vertices at zeta = -1 + * - Node 3-5: triangle vertices at zeta = +1 (node a + 3 is directly above node a) + * - Nodes 6 + 2a and 7+2a , a = 0, 1, 2 : on the vertical edge through vertex a, + * at zeta = -1/3 and zeta = +1/3 + * + * Barycentric coordinates of the triangle: L0 = 1 - xi0 - xi1, L1 = xi0, L2 = xi1. + * Cubic Lagrange functions in zeta on the stations -1, -1/3, 1/3, 1: + * \f[ + * \ell_0 = -\tfrac{9}{16}(\zeta+\tfrac{1}{3})(\zeta-\tfrac{1}{3})(\zeta-1), \qquad + * \ell_1 = \tfrac{27}{16}(\zeta+1)(\zeta-\tfrac{1}{3})(\zeta-1), + * \f] + * \f[ + * \ell_2 = -\tfrac{27}{16}(\zeta+1)(\zeta+\tfrac{1}{3})(\zeta-1), \qquad + * \ell_3 = \tfrac{9}{16}(\zeta+1)(\zeta+\tfrac{1}{3})(\zeta-\tfrac{1}{3}). + * \f] + * Shape functions, for a = 0, 1, 2: + * - \f$N_a = L_a \, \ell_0\f$ and \f$N_{a+3} = L_a \, \ell_3\f$ + * - \f$N_{6+2a} = L_a \, \ell_1\f$ and \f$N_{7+2a} = L_a \, \ell_2\f$ + * + */ + +struct LinearCubicWedgeShape { + // clang-format off + static const size_t numNodes = 12; ///< Number of nodes (12 for linear-cubic wedge) + static const size_t meshEntDim = 3; ///< Mesh entity dimension (3D) + constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge}; ///< DOFs at vertices and edges + constexpr static size_t NumDofHolders[2] = {6, 3}; ///< 6 vertices, 3 edges + constexpr static size_t DofsPerHolder[2] = {1, 2}; ///< 1 DOF per vertex, 2 DOF per edge + constexpr static size_t Order = 3; ///< Polynomial order (cubic) + // 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 { + // nodes at vertices + 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 + + // nodes at two stations of edges at -1/3 and +1/3 + 0, 0, -1.0/3, // node 6 + 0, 0, 1.0/3, // node 7 + 1, 0, -1.0/3, // node 8 + 1, 0, 1.0/3, // node 9 + 0, 1, -1.0/3, // node 10 + 0, 1, 1.0/3 // node 11 + }; + // clang-format on + } + + /** + * @brief Evaluate shape functions at parametric coordinate + * @param xi Parametric coordinates (xi0, xi1, zeta) where xi0, xi1 + * \f$\ge\f$ 0, xi0 + xi1 \f$\le\f$ 1 and zeta \f$\in\f$ [-1,1] + * @return Array of shape function values [N0, ..., N11] where, for a = 0, 1, 2, + * \f$N_a = L_a \, \ell_0(\zeta)\f$, \f$N_{a+3} = L_a \, \ell_3(\zeta)\f$, + * \f$N_{6+2a} = L_a \, \ell_1(\zeta)\f$ and \f$N_{7+2a} = L_a \, \ell_2(\zeta)\f$, + * with L0 = 1 - xi0 - xi1, L1 = xi0, L2 = xi1 and the cubic Lagrange + * functions \f$\ell_0, \dots, \ell_3\f$ defined in the class description + */ + 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)); + assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); + const Real L1 = xi[0]; + const Real L2 = xi[1]; + + const Real l0 = -1.0/16 * (3 * xi[2] + 1) * (3 * xi[2] - 1) * (xi[2] - 1); //cubic edge function at -1 station + const Real l1 = 9.0/16 * (xi[2] + 1) * (3 * xi[2] - 1) * (xi[2] - 1); // cubic edge function at -1/3 station + const Real l2 = -9.0/16 * (xi[2] + 1) * (3 * xi[2] + 1) * (xi[2] - 1); // cubic edge function at +1/3 station + const Real l3 = 1.0/16 * (xi[2] + 1) * (3 * xi[2] + 1) * (3 * xi[2] - 1); // cubic edge function at +1 station + + // clang-format off + return {L0 * l0, + L1 * l0, + L2 * l0, + L0 * l3, + L1 * l3, + L2 * l3, + L0 * l1, + L0 * l2, + L1 * l1, + L1 * l2, + L2 * l1, + L2 * l2 + }; + // 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 l0 = -1.0/16 * (3 * xi[2] + 1) * (3 * xi[2] - 1) * (xi[2] - 1); //cubic edge function at -1 station + const Real dl0 = 1.0/16 * (-27 * xi[2] * xi[2] + 18 * xi[2] + 1); //derivative of cubic edge function at -1 station + + const Real l1 = 9.0/16 * (xi[2] + 1) * (3 * xi[2] - 1) * (xi[2] - 1); // cubic edge function at -1/3 station + const Real dl1 = 1.0/16 * (81 * xi[2] * xi[2] - 18 * xi[2] -27);// derivative of cubic edge function at -1/3 station + + const Real l2 = -9.0/16 * (xi[2] + 1) * (3 * xi[2] + 1) * (xi[2] - 1); // cubic edge function at +1/3 station + const Real dl2 = 1.0/16 * (-81 * xi[2] * xi[2] - 18 * xi[2] + 27); // derivative of cubic edge function at +1/3 station + + const Real l3 = 1.0/16 * (xi[2] + 1) * (3 * xi[2] + 1) * (3 * xi[2] - 1); // cubic edge function at +1 station + const Real dl3 = 1.0/16 * (27 * xi[2] * xi[2] + 18 * xi[2] - 1); // derivative of cubic edge function at +1 station + + + // clang-format off + return {-l0, -l0, L0 * dl0, //node 0 + l0, 0, L1 * dl0, //node 1 + 0, l0, L2 * dl0, //node 2 + -l3, -l3, L0 * dl3, //node 3 + l3, 0, L1 * dl3, //node 4 + 0, l3, L2 * dl3, //node 5 + // nodes at two stations of edges at -1/3 and +1/3 + -l1, -l1, L0 * dl1, // node 6 + -l2, -l2, L0 * dl2, // node 7 + l1, 0 , L1 * dl1, // node 8 + l2, 0 , L1 * dl2, // node 9 + 0, l1, L2 * dl1, // node 10 + 0, l2 , L2 * dl2 // node 11 + }; + // clang-format on + } +}; + } // namespace MeshField #endif From fd39ae6d5baa739b88c0204895e663b24b5325ee Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 15 Oct 2026 13:56:13 -0400 Subject: [PATCH 11/23] put all the tests inside a single function --- test/testCubicWedgeElement.cpp | 213 +++++++++++++++++++++++++++++++++ test/testCubicWedgeShape.cpp | 148 +++++++++++++++++++++++ 2 files changed, 361 insertions(+) create mode 100644 test/testCubicWedgeElement.cpp create mode 100644 test/testCubicWedgeShape.cpp diff --git a/test/testCubicWedgeElement.cpp b/test/testCubicWedgeElement.cpp new file mode 100644 index 0000000..3a968c5 --- /dev/null +++ b/test/testCubicWedgeElement.cpp @@ -0,0 +1,213 @@ +#include +#include +#include "MeshField_Defines.hpp" +#include "MeshField_Shape.hpp" +#include "MeshField_Element.hpp" +#include +#include +#include +#include + +// Checks MeshField::LinearCubicWedgeShape inside a MeshField::FieldElement: +// field evaluation (MeshField::evaluate) and the Jacobian +// (MeshField::getJacobians, MeshField::getJacobianDeterminants) for two cubic +// wedges. The nodes are placed by geometric maps that are linear in (xi0, xi1) +// and cubic in zeta, so the element map reproduces them exactly and the +// expected values are the analytic maps and their derivatives. Failures are +// counted explicitly so the test also fails in builds where assert() is +// disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; +using Vec3 = Omega_h::Vector<3>; + +namespace { + +constexpr int numWedges = 2; +constexpr int numNodes = MeshField::LinearCubicWedgeShape::numNodes; +constexpr int numVtxPerWedge = 6; +constexpr int numEdgesPerWedge = 3; + +bool isClose(Real a, Real b) { return Omega_h::are_close(a, b, 1e-12, 1e-12); } + +// Field storage for the test: one value per vertex and two per vertical edge, +// each with three components. FieldElement only calls +// field(entity, dof, comp, topo), so this stands in for a MeshField field. +struct CubicWedgeField { + using BaseType = Real; + static constexpr size_t numComp = 3; + Kokkos::View vtx; // (vertex, comp) + Kokkos::View edge; // (edge, dof, comp) + KOKKOS_INLINE_FUNCTION Real &operator()(MeshField::LO ent, MeshField::LO dof, + MeshField::LO comp, + MeshField::Mesh_Topology topo) const { + return topo == MeshField::Vertex ? vtx(ent, comp) : edge(ent, dof, comp); + } +}; + +// Local node n of wedge w: nodes 0-5 are vertex 6*w + n; nodes 6+2a and 7+2a +// are the two DOFs on vertical edge 3*w + a, at zeta = -1/3 and +1/3. +struct CubicWedgeToField { + 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); + if (wedgeNodeIdx < numVtxPerWedge) { + const MeshField::LO vtx = numVtxPerWedge * wedge + wedgeNodeIdx; + return {0, wedgeCompIdx, vtx, MeshField::Vertex}; + } + const MeshField::LO a = (wedgeNodeIdx - numVtxPerWedge) / 2; + const MeshField::LO dof = (wedgeNodeIdx - numVtxPerWedge) % 2; + const MeshField::LO edge = numEdgesPerWedge * wedge + a; + return {dof, wedgeCompIdx, edge, MeshField::Edge}; + } +}; + +// Geometric maps that place the nodes of each wedge. Both are linear in +// (xi0, xi1) and cubic in zeta, so they lie in the element's function space. +Vec3 mapX(int w, Vec3 xi) { + const Real x0 = xi[0], x1 = xi[1], z = xi[2]; + if (w == 0) + return {2 * x0 + 0.3 * z * z * x1 + 0.2 * z * z * z, + 3 * x1 + 0.1 * z * z * z * x0 - 0.4 * z, + 2 + 2 * z + 0.25 * z * z * z}; + return {x0 + 0.5 * x1 + 0.3 * z * z * z * x0 + 0.1 * z, + x1 + 0.2 * z * z * x1 + 0.4 * z * z, + 1.5 + z + 0.2 * z * z * z * x1 + 0.1 * z * z * z}; +} + +// Analytic derivatives of mapX. Row i is dx/dxi_i, the same layout as the +// MeshField Jacobian: J(i, j) = dx_j / dxi_i. +std::array mapJ(int w, Vec3 xi) { + const Real x0 = xi[0], x1 = xi[1], z = xi[2]; + if (w == 0) + return {Vec3{2, 0.1 * z * z * z, 0}, + Vec3{0.3 * z * z, 3, 0}, + Vec3{0.6 * z * x1 + 0.6 * z * z, 0.3 * z * z * x0 - 0.4, + 2 + 0.75 * z * z}}; + return {Vec3{1 + 0.3 * z * z * z, 0, 0}, + Vec3{0.5, 1 + 0.2 * z * z, 0.2 * z * z * z}, + Vec3{0.9 * z * z * x0 + 0.1, 0.4 * z * x1 + 0.8 * z, + 1 + 0.6 * z * z * x1 + 0.3 * z * z}}; +} + +// Reference coordinates of the 12 nodes: six vertices, then two edge nodes on +// each vertical edge at zeta = -1/3 and +1/3. +const Vec3 nodeRef[numNodes] = { + {0, 0, -1}, {1, 0, -1}, {0, 1, -1}, {0, 0, 1}, + {1, 0, 1}, {0, 1, 1}, {0, 0, -1. / 3}, {0, 0, 1. / 3}, + {1, 0, -1. / 3}, {1, 0, 1. / 3}, {0, 1, -1. / 3}, {0, 1, 1. / 3}}; + +// The 12 nodes, then interior points. +std::vector makeLocalPoints() { + std::vector pts(nodeRef, nodeRef + numNodes); + pts.push_back({1. / 3, 1. / 3, 0}); + pts.push_back({0.2, 0.3, -0.4}); + pts.push_back({0.6, 0.1, 0.7}); + pts.push_back({0.05, 0.9, 0.95}); + return pts; +} + +// Vertex and edge DOFs, filled from the geometric maps at the nodes. The +// device views are filled through their mirrors, so the layouts always match. +CubicWedgeField makeCubicWedgeField() { + const int numVtx = numWedges * numVtxPerWedge; + const int numEdges = numWedges * numEdgesPerWedge; + Kokkos::View vtx("vtx", numVtx, 3); + Kokkos::View edge("edge", numEdges, 2, 3); + auto vtx_h = Kokkos::create_mirror_view(vtx); + auto edge_h = Kokkos::create_mirror_view(edge); + for (int w = 0; w < numWedges; ++w) + for (int n = 0; n < numNodes; ++n) { + const Vec3 X = mapX(w, nodeRef[n]); + for (int d = 0; d < 3; ++d) { + if (n < numVtxPerWedge) { + vtx_h(numVtxPerWedge * w + n, d) = X[d]; + } else { + const int a = (n - numVtxPerWedge) / 2; + const int dof = (n - numVtxPerWedge) % 2; + edge_h(numEdgesPerWedge * w + a, dof, d) = X[d]; + } + } + } + Kokkos::deep_copy(vtx, vtx_h); + Kokkos::deep_copy(edge, edge_h); + return CubicWedgeField{vtx, edge}; +} + +// Compares the evaluated positions, Jacobians and determinants of wedge w +// with the analytic maps at every point in pts. Returns the failure count. +template +int checkWedge(int w, const std::vector &pts, const XView &x_h, + const JView &J_h, const DetView &detJ_h) { + const int numPts = static_cast(pts.size()); + int evalFailed = 0, jacFailed = 0, detFailed = 0; + for (int p = 0; p < numPts; ++p) { + const int row = w * numPts + p; + const Vec3 xi = pts[p]; + const Vec3 xExpected = mapX(w, xi); + for (int d = 0; d < 3; ++d) + evalFailed += !isClose(x_h(row, d), xExpected[d]); + const auto jExpected = mapJ(w, xi); + 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 << "LinearCubicWedgeShape FieldElement, cubic wedge " << w + << ": evaluate failures " << evalFailed << ", Jacobian failures " + << jacFailed << ", det J failures " << detFailed << "\n"; + return evalFailed + jacFailed + detFailed; +} + +// Builds the field and the FieldElement, runs the evaluation, Jacobian and +// determinant, and checks every wedge. +int runCubicWedgeElementChecks() { + const auto pts = makeLocalPoints(); + const int numPts = static_cast(pts.size()); + + CubicWedgeField field = makeCubicWedgeField(); + MeshField::FieldElement fe(numWedges, field, + MeshField::LinearCubicWedgeShape(), + CubicWedgeToField()); + + // 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) = pts[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, pts, x_h, J_h, detJ_h); + return numFailed; +} + +} // namespace + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + const int numFailed = runCubicWedgeElementChecks(); + Kokkos::finalize(); + return numFailed ? 1 : 0; +} diff --git a/test/testCubicWedgeShape.cpp b/test/testCubicWedgeShape.cpp new file mode 100644 index 0000000..0a3207e --- /dev/null +++ b/test/testCubicWedgeShape.cpp @@ -0,0 +1,148 @@ +#include +#include +#include +#include + +// Checks of MeshField::LinearCubicWedgeShape, linear in (xi0, xi1) and cubic +// in zeta. 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. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; +using Shape = MeshField::LinearCubicWedgeShape; + +KOKKOS_INLINE_FUNCTION bool isClose(Real a, Real b, Real tol) { + return Kokkos::fabs(a - b) <= tol; +} + +// A function in the span of the cubic wedge basis, +// {1, xi0, xi1} x {1, zeta, zeta^2, zeta^3}, and its gradient +KOKKOS_INLINE_FUNCTION Real f(Real x0, Real x1, Real z) { + return 0.5 + 2 * x0 - 3 * x1 + 0.25 * z - 1.5 * z * z + 0.75 * z * z * z + + 4 * x0 * z - 1.5 * x1 * z * z - 2 * x0 * z * z * z + x1 * z * z * z; +} +KOKKOS_INLINE_FUNCTION MeshField::Vector3 gradF(Real x0, Real x1, Real z) { + return {2 + 4 * z - 2 * z * z * z, + -3 - 1.5 * z * z + z * z * z, + 0.25 - 3 * z + 2.25 * z * z + 4 * x0 - 3 * x1 * z - 6 * x0 * z * z + + 3 * x1 * z * z}; +} + +// Each shape function is one at its own node and zero at the others. +int checkKroneckerDelta() { + int numFailed = 0; + Kokkos::parallel_reduce( + "cubicWedgeKronecker", + 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, 10 * MeshField::MachinePrecision); + } + }, + numFailed); + return numFailed; +} + +// Interior points, including points at the zeta stations and between them. +Kokkos::View makeInteriorPoints() { + // 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}, + {0.2, 0.3, 1. / 3}, // at a station + {0.6, 0.1, -1. / 3}, // at a station + {0.3, 0.2, -0.7}, // between -1 and -1/3 + {0.1, 0.8, 0.5}}; // between 1/3 and 1 + // clang-format on + const int numPts = static_cast(sizeof(pts) / sizeof(pts[0])); + Kokkos::View pts_d("cubicWedgePts", 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, 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( + "cubicWedgeInterior", Kokkos::RangePolicy(0, numPts), + KOKKOS_LAMBDA(const int i, int &failed) { + const Real tol = 10 * MeshField::MachinePrecision; + 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; + } + } + 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); + } + + // gradients against central finite differences of the values + 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; +} + +// Runs both checks and prints their failure counts. +int runCubicWedgeShapeChecks() { + const int kronecker = checkKroneckerDelta(); + std::cout << "LinearCubicWedgeShape Kronecker delta failures: " << kronecker + << "\n"; + const int interior = checkInteriorPoints(); + std::cout << "LinearCubicWedgeShape interior point failures: " << interior + << "\n"; + return kronecker + interior; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + const int numFailed = runCubicWedgeShapeChecks(); + Kokkos::finalize(); + return numFailed ? 1 : 0; +} From 230d2d8707f32d0592b9b09ab0e6090b55cae419 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 15 Oct 2026 13:56:48 -0400 Subject: [PATCH 12/23] add testCubicWedgeShape.cpp and testCubicWedgeElement.cpp test files --- CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index b60c19c..f2cc527 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -214,7 +214,9 @@ 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(CubicWedgeShapeTest test/testCubicWedgeShape.cpp) meshfields_add_exe(WedgeElementTest test/testWedgeElement.cpp) +meshfields_add_exe(CubicWedgeElementTest test/testCubicWedgeElement.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) if(MeshFields_USE_Cabana) @@ -238,7 +240,9 @@ smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(OmegahTetTest ./OmegahTetTest) test_func(IntegratorMeasure ./IntegratorMeasure) test_func(WedgeShapeTest ./WedgeShapeTest) +test_func(CubicWedgeShapeTest ./CubicWedgeShapeTest) test_func(WedgeElementTest ./WedgeElementTest) +test_func(CubicWedgeElementTest ./CubicWedgeElementTest) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) 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 13/23] 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 14/23] 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 b67403e1f76d89104b02a3da8e4e7ffcd71ce092 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 15 Oct 2026 14:47:06 -0400 Subject: [PATCH 15/23] remove arrows --- CMakeLists.txt | 1 - test/testLinearWedgeElement.cpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 35f7a97..44f6685 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -217,7 +217,6 @@ meshfields_add_exe(CubicWedgeShapeTest test/testCubicWedgeShape.cpp) meshfields_add_exe(CubicWedgeElementTest test/testCubicWedgeElement.cpp) meshfields_add_exe(LinearWedgeShapeTest test/testLinearWedgeShape.cpp) meshfields_add_exe(LinearWedgeElementTest test/testLinearWedgeElement.cpp) ->>>>>>> linear-wedge meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) if(MeshFields_USE_Cabana) diff --git a/test/testLinearWedgeElement.cpp b/test/testLinearWedgeElement.cpp index 0d46d3f..7b2638f 100644 --- a/test/testLinearWedgeElement.cpp +++ b/test/testLinearWedgeElement.cpp @@ -190,4 +190,4 @@ int main(int argc, char **argv) { const int numFailed = runLinearWedgeElementChecks(); Kokkos::finalize(); return numFailed ? 1 : 0; -} \ No newline at end of file +} From 8c596a6635f1e73aed9ea943ed3b38da45239462 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 13:42:53 -0400 Subject: [PATCH 16/23] rename from testCubicWedgeShape.cpp to testLinearCubicWedgeShape.cpp --- test/{testCubicWedgeShape.cpp => testLinearCubicWedgeShape.cpp} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/{testCubicWedgeShape.cpp => testLinearCubicWedgeShape.cpp} (100%) diff --git a/test/testCubicWedgeShape.cpp b/test/testLinearCubicWedgeShape.cpp similarity index 100% rename from test/testCubicWedgeShape.cpp rename to test/testLinearCubicWedgeShape.cpp From b67a673958282b1f95e3697771c6d19861a90585 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 13:43:43 -0400 Subject: [PATCH 17/23] rename testCubicWedgeElement.cpp to testLinearCubicWedgeElement.cpp --- ...{testCubicWedgeElement.cpp => testLinearCubicWedgeElement.cpp} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/{testCubicWedgeElement.cpp => testLinearCubicWedgeElement.cpp} (100%) diff --git a/test/testCubicWedgeElement.cpp b/test/testLinearCubicWedgeElement.cpp similarity index 100% rename from test/testCubicWedgeElement.cpp rename to test/testLinearCubicWedgeElement.cpp From aa13b618d1373849d9e59852304bb952a66da82f Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 13:45:37 -0400 Subject: [PATCH 18/23] rename CubicWedge to LinearCubicWedge --- CMakeLists.txt | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 44f6685..0771ae3 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(CubicWedgeShapeTest test/testCubicWedgeShape.cpp) -meshfields_add_exe(CubicWedgeElementTest test/testCubicWedgeElement.cpp) +meshfields_add_exe(LinearCubicWedgeShapeTest test/testLinearCubicWedgeShape.cpp) +meshfields_add_exe(LinearCubicWedgeElementTest test/testLinearCubicWedgeElement.cpp) meshfields_add_exe(LinearWedgeShapeTest test/testLinearWedgeShape.cpp) meshfields_add_exe(LinearWedgeElementTest test/testLinearWedgeElement.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) @@ -239,8 +239,8 @@ test_func(CountIntegrator ./CountIntegrator) smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(OmegahTetTest ./OmegahTetTest) test_func(IntegratorMeasure ./IntegratorMeasure) -test_func(CubicWedgeShapeTest ./CubicWedgeShapeTest) -test_func(CubicWedgeElementTest ./CubicWedgeElementTest) +test_func(LinearCubicWedgeShapeTest ./LinearCubicWedgeShapeTest) +test_func(LinearCubicWedgeElementTest ./LinearCubicWedgeElementTest) test_func(LinearWedgeShapeTest ./LinearWedgeShapeTest) test_func(LinearWedgeElementTest ./LinearWedgeElementTest) if(MeshFields_USE_EXCEPTIONS) From 649334886780d743b7a3f14f29ff0873281738c1 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Mon, 26 Oct 2026 16:56:12 -0400 Subject: [PATCH 19/23] Add LinearCubicWedgeShape storage to CreateLagrangeField --- src/MeshField_ShapeField.hpp | 88 ++++++++++++++++++++++++++++++++---- 1 file changed, 80 insertions(+), 8 deletions(-) diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 90ab0ea..9133400 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -8,7 +8,7 @@ #include "MeshField_Field.hpp" #include "MeshField_Shape.hpp" #include //decltype -#include // std::forward +#include // std::forward namespace MeshField { @@ -83,8 +83,8 @@ struct ShapeField : public Mixins... { // sites like ShapeField(meshInfo, {vtxField}) are fine; passing a named // lvalue requires std::move). mixins is moved-from by std::forward below // - do not read mixins again after the initializer list. - ShapeField(const MeshInfo &meshInfoIn, Mixins &&... mixins) - : meshInfo(meshInfoIn), Mixins( std::forward(mixins) )... {}; + ShapeField(const MeshInfo &meshInfoIn, Mixins &&...mixins) + : meshInfo(meshInfoIn), Mixins(std::forward(mixins))... {}; }; //! [ShapeField] @@ -157,6 +157,17 @@ template struct LinearAccessor { } }; //! [LinearAccessor] +/** + * @brief + * Shape functions used by CreateLagrangeField when the caller does not name + * one: triangles for 1d and 2d meshes and tetrahedra for 3d meshes. + */ +template +using DefaultLagrangeShape = std::conditional_t< + order == 2, + std::conditional_t, + std::conditional_t>; /** * @brief @@ -174,6 +185,9 @@ template struct LinearAccessor { * @tparam order the order of the shape functions; linear(1) and quadratic(2) * are supported * @tparam dim the dimension of the mesh + * @tparam Shape the shape functions; the default picks the triangle or + * tetrahedron shape, and LinearCubicWedgeShape (order 3) stores one DOF per + * vertex and two per vertical edge * * @param meshInfo defines on-process mesh metadata * @return a linear or quadratic ShapeField @@ -181,19 +195,75 @@ template struct LinearAccessor { template typename Controller = MeshField::KokkosController, - typename DataType, size_t order, size_t dim, size_t numComp> + typename DataType, size_t order, size_t dim, size_t numComp, + typename Shape = DefaultLagrangeShape> auto CreateLagrangeField(const MeshInfo &meshInfo) { + constexpr bool isLinearCubicWedge = + std::is_same_v; static_assert((std::is_same_v == true || std::is_same_v == true), "CreateLagrangeField only supports single and double precision " "floating point fields\n"); - static_assert( - (order == 1 || order == 2), - "CreateLagrangeField only supports linear and quadratic fields\n"); + static_assert((order == 1 || order == 2 || isLinearCubicWedge), + "CreateLagrangeField only supports linear and quadratic " + "fields, and LinearCubicWedgeShape\n"); + static_assert(isLinearCubicWedge || + std::is_same_v>, + "CreateLagrangeField: Shape must be the default or " + "LinearCubicWedgeShape\n"); + static_assert(Shape::Order == order, + "CreateLagrangeField order must match Shape::Order\n"); static_assert((dim == 1 || dim == 2 || dim == 3), "CreateLagrangeField only supports 1d, 2d, and 3d meshes\n"); using MemorySpace = typename ExecutionSpace::memory_space; - if constexpr (order == 1 && (dim == 1 || dim == 2 || dim == 3)) { + if constexpr (isLinearCubicWedge) { + static_assert(dim == 3, "LinearCubicWedgeShape fields need a 3d mesh\n"); + if (meshInfo.numVtx <= 0) { + fail("mesh has no vertices\n"); + } + if (meshInfo.numEdge <= 0) { + fail("mesh has no edges\n"); + } + if (meshInfo.numWedge <= 0) { + fail("mesh has no wedges\n"); + } + // 1 dof per vertex and Shape::DofsPerHolder[1] (2) dofs per vertical edge + constexpr int dofsPerEdge = Shape::DofsPerHolder[1]; +#ifdef MESHFIELDS_ENABLE_CABANA + using Ctrlr = std::conditional_t< + std::is_same_v< + Controller, + MeshField::CabanaController>, + Controller, + Controller>; + auto createController = [](auto numVtx, auto numEdge) { + if constexpr (std::is_same_v< + Controller, + MeshField::CabanaController>) { + return Ctrlr({numVtx, numEdge}); + } else { + return Ctrlr({/*field 0*/ numVtx, 1, numComp, + /*field 1*/ numEdge, dofsPerEdge, numComp}); + } + }; + Ctrlr kk_ctrl = createController(meshInfo.numVtx, meshInfo.numEdge); +#else + using Ctrlr = + Controller; + Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp, + /*field 1*/ meshInfo.numEdge, dofsPerEdge, numComp}); +#endif + auto vtxField = MeshField::makeField(kk_ctrl); + auto edgeField = MeshField::makeField(kk_ctrl); + // QuadraticAccessor forwards the dof index to the edge field, so it also + // serves edges with more than one dof + using QA = QuadraticAccessor; + using WedgeShapeField = ShapeField; + WedgeShapeField wsf(meshInfo, {vtxField, edgeField}); + return FieldWithController{kk_ctrl, wsf}; + } else if constexpr (order == 1 && (dim == 1 || dim == 2 || dim == 3)) { if (meshInfo.numVtx <= 0) { fail("mesh has no vertices\n"); } @@ -315,3 +385,5 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { } // namespace MeshField #endif + + From dae3656276880a47edca0f7bdddf1168ecf2c50e Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Mon, 26 Oct 2026 17:12:29 -0400 Subject: [PATCH 20/23] Store LinearCubicWedgeElementTest coordinates in a library field --- test/testLinearCubicWedgeElement.cpp | 97 ++++++++++++++-------------- 1 file changed, 50 insertions(+), 47 deletions(-) diff --git a/test/testLinearCubicWedgeElement.cpp b/test/testLinearCubicWedgeElement.cpp index 3a968c5..78b9242 100644 --- a/test/testLinearCubicWedgeElement.cpp +++ b/test/testLinearCubicWedgeElement.cpp @@ -1,8 +1,9 @@ -#include -#include #include "MeshField_Defines.hpp" -#include "MeshField_Shape.hpp" #include "MeshField_Element.hpp" +#include "MeshField_Shape.hpp" +#include "MeshField_ShapeField.hpp" +#include +#include #include #include #include @@ -11,11 +12,10 @@ // Checks MeshField::LinearCubicWedgeShape inside a MeshField::FieldElement: // field evaluation (MeshField::evaluate) and the Jacobian // (MeshField::getJacobians, MeshField::getJacobianDeterminants) for two cubic -// wedges. The nodes are placed by geometric maps that are linear in (xi0, xi1) -// and cubic in zeta, so the element map reproduces them exactly and the -// expected values are the analytic maps and their derivatives. Failures are -// counted explicitly so the test also fails in builds where assert() is -// disabled. +// wedges. The coordinates are stored in a field from +// MeshField::CreateLagrangeField with LinearCubicWedgeShape: one DOF per vertex +// and two per vertical edge. +// The nodes are placed by geometric maps that are linear in (xi0, xi1) using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; @@ -30,21 +30,6 @@ constexpr int numEdgesPerWedge = 3; bool isClose(Real a, Real b) { return Omega_h::are_close(a, b, 1e-12, 1e-12); } -// Field storage for the test: one value per vertex and two per vertical edge, -// each with three components. FieldElement only calls -// field(entity, dof, comp, topo), so this stands in for a MeshField field. -struct CubicWedgeField { - using BaseType = Real; - static constexpr size_t numComp = 3; - Kokkos::View vtx; // (vertex, comp) - Kokkos::View edge; // (edge, dof, comp) - KOKKOS_INLINE_FUNCTION Real &operator()(MeshField::LO ent, MeshField::LO dof, - MeshField::LO comp, - MeshField::Mesh_Topology topo) const { - return topo == MeshField::Vertex ? vtx(ent, comp) : edge(ent, dof, comp); - } -}; - // Local node n of wedge w: nodes 0-5 are vertex 6*w + n; nodes 6+2a and 7+2a // are the two DOFs on vertical edge 3*w + a, at zeta = -1/3 and +1/3. struct CubicWedgeToField { @@ -86,8 +71,7 @@ Vec3 mapX(int w, Vec3 xi) { std::array mapJ(int w, Vec3 xi) { const Real x0 = xi[0], x1 = xi[1], z = xi[2]; if (w == 0) - return {Vec3{2, 0.1 * z * z * z, 0}, - Vec3{0.3 * z * z, 3, 0}, + return {Vec3{2, 0.1 * z * z * z, 0}, Vec3{0.3 * z * z, 3, 0}, Vec3{0.6 * z * x1 + 0.6 * z * z, 0.3 * z * z * x0 - 0.4, 2 + 0.75 * z * z}}; return {Vec3{1 + 0.3 * z * z * z, 0, 0}, @@ -99,9 +83,9 @@ std::array mapJ(int w, Vec3 xi) { // Reference coordinates of the 12 nodes: six vertices, then two edge nodes on // each vertical edge at zeta = -1/3 and +1/3. const Vec3 nodeRef[numNodes] = { - {0, 0, -1}, {1, 0, -1}, {0, 1, -1}, {0, 0, 1}, - {1, 0, 1}, {0, 1, 1}, {0, 0, -1. / 3}, {0, 0, 1. / 3}, - {1, 0, -1. / 3}, {1, 0, 1. / 3}, {0, 1, -1. / 3}, {0, 1, 1. / 3}}; + {0, 0, -1}, {1, 0, -1}, {0, 1, -1}, {0, 0, 1}, + {1, 0, 1}, {0, 1, 1}, {0, 0, -1. / 3}, {0, 0, 1. / 3}, + {1, 0, -1. / 3}, {1, 0, 1. / 3}, {0, 1, -1. / 3}, {0, 1, 1. / 3}}; // The 12 nodes, then interior points. std::vector makeLocalPoints() { @@ -113,15 +97,15 @@ std::vector makeLocalPoints() { return pts; } -// Vertex and edge DOFs, filled from the geometric maps at the nodes. The -// device views are filled through their mirrors, so the layouts always match. -CubicWedgeField makeCubicWedgeField() { +// Writes the node positions from mapX into the coordinate field: vertex +// nodes to the vertex DOFs, and the zeta = -1/3 and +1/3 nodes to DOFs 0 and 1 +// of their vertical edge. +template void setCubicWedgeCoords(CoordField field) { const int numVtx = numWedges * numVtxPerWedge; const int numEdges = numWedges * numEdgesPerWedge; - Kokkos::View vtx("vtx", numVtx, 3); - Kokkos::View edge("edge", numEdges, 2, 3); - auto vtx_h = Kokkos::create_mirror_view(vtx); - auto edge_h = Kokkos::create_mirror_view(edge); + Kokkos::View vtx_h("vtx_h", numVtx); + Kokkos::View edge_h("edge_h", numEdges); + for (int w = 0; w < numWedges; ++w) for (int n = 0; n < numNodes; ++n) { const Vec3 X = mapX(w, nodeRef[n]); @@ -135,9 +119,22 @@ CubicWedgeField makeCubicWedgeField() { } } } - Kokkos::deep_copy(vtx, vtx_h); - Kokkos::deep_copy(edge, edge_h); - return CubicWedgeField{vtx, edge}; + + auto vtx = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), vtx_h); + auto edge = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), edge_h); + Kokkos::parallel_for( + "setVtxCoords", Kokkos::RangePolicy(0, numVtx), + KOKKOS_LAMBDA(const int v) { + for (int d = 0; d < 3; ++d) + field(v, 0, d, MeshField::Vertex) = vtx(v, d); + }); + Kokkos::parallel_for( + "setEdgeCoords", Kokkos::RangePolicy(0, numEdges), + KOKKOS_LAMBDA(const int e) { + for (int dof = 0; dof < 2; ++dof) + for (int d = 0; d < 3; ++d) + field(e, dof, d, MeshField::Edge) = edge(e, dof, d); + }); } // Compares the evaluated positions, Jacobians and determinants of wedge w @@ -168,16 +165,24 @@ int checkWedge(int w, const std::vector &pts, const XView &x_h, } // Builds the field and the FieldElement, runs the evaluation, Jacobian and -// determinant, and checks every wedge. +// determinant, and checks every wedge. int runCubicWedgeElementChecks() { const auto pts = makeLocalPoints(); const int numPts = static_cast(pts.size()); - - CubicWedgeField field = makeCubicWedgeField(); - MeshField::FieldElement fe(numWedges, field, + // numEdge counts only the vertical edges, which hold the zeta = +-1/3 DOFs + const MeshField::MeshInfo meshInfo{.numVtx = numWedges * numVtxPerWedge, + .numEdge = numWedges * numEdgesPerWedge, + .numWedge = numWedges, + .dim = 3}; + auto coordField = + MeshField::CreateLagrangeField( + meshInfo); + setCubicWedgeCoords(coordField.field); + MeshField::FieldElement fe(numWedges, coordField.field, MeshField::LinearCubicWedgeShape(), CubicWedgeToField()); - // 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); @@ -190,10 +195,8 @@ int runCubicWedgeElementChecks() { 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 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); From 5cf25f847ff0238ae7dc69858ecfe668e5f62948 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 02:29:58 -0400 Subject: [PATCH 21/23] Restore test comment and apply clang-format to the cubic shape test --- src/MeshField_ShapeField.hpp | 9 ++++----- test/testLinearCubicWedgeElement.cpp | 4 ++++ test/testLinearCubicWedgeShape.cpp | 5 ++--- test/testLinearWedgeElement.cpp | 2 +- 4 files changed, 11 insertions(+), 9 deletions(-) diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 9133400..cfb4502 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -8,7 +8,7 @@ #include "MeshField_Field.hpp" #include "MeshField_Shape.hpp" #include //decltype -#include // std::forward +#include // std::forward namespace MeshField { @@ -83,8 +83,8 @@ struct ShapeField : public Mixins... { // sites like ShapeField(meshInfo, {vtxField}) are fine; passing a named // lvalue requires std::move). mixins is moved-from by std::forward below // - do not read mixins again after the initializer list. - ShapeField(const MeshInfo &meshInfoIn, Mixins &&...mixins) - : meshInfo(meshInfoIn), Mixins(std::forward(mixins))... {}; + ShapeField(const MeshInfo &meshInfoIn, Mixins &&... mixins) + : meshInfo(meshInfoIn), Mixins( std::forward(mixins) )... {}; }; //! [ShapeField] @@ -157,6 +157,7 @@ template struct LinearAccessor { } }; //! [LinearAccessor] + /** * @brief * Shape functions used by CreateLagrangeField when the caller does not name @@ -385,5 +386,3 @@ auto CreateCoordinateField(const MeshInfo &meshInfo) { } // namespace MeshField #endif - - diff --git a/test/testLinearCubicWedgeElement.cpp b/test/testLinearCubicWedgeElement.cpp index 78b9242..0ec2c17 100644 --- a/test/testLinearCubicWedgeElement.cpp +++ b/test/testLinearCubicWedgeElement.cpp @@ -16,6 +16,10 @@ // MeshField::CreateLagrangeField with LinearCubicWedgeShape: one DOF per vertex // and two per vertical edge. // The nodes are placed by geometric maps that are linear in (xi0, xi1) +// and cubic in zeta, so the element map reproduces them exactly and the +// expected values are the analytic maps and their derivatives. Failures are +// counted explicitly so the test also fails in builds where assert() is +// disabled. using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; diff --git a/test/testLinearCubicWedgeShape.cpp b/test/testLinearCubicWedgeShape.cpp index 0a3207e..8b37471 100644 --- a/test/testLinearCubicWedgeShape.cpp +++ b/test/testLinearCubicWedgeShape.cpp @@ -23,8 +23,7 @@ KOKKOS_INLINE_FUNCTION Real f(Real x0, Real x1, Real z) { 4 * x0 * z - 1.5 * x1 * z * z - 2 * x0 * z * z * z + x1 * z * z * z; } KOKKOS_INLINE_FUNCTION MeshField::Vector3 gradF(Real x0, Real x1, Real z) { - return {2 + 4 * z - 2 * z * z * z, - -3 - 1.5 * z * z + z * z * z, + return {2 + 4 * z - 2 * z * z * z, -3 - 1.5 * z * z + z * z * z, 0.25 - 3 * z + 2.25 * z * z + 4 * x0 - 3 * x1 * z - 6 * x0 * z * z + 3 * x1 * z * z}; } @@ -129,7 +128,7 @@ int checkInteriorPoints() { return numFailed; } -// Runs both checks and prints their failure counts. +// Runs both checks and prints their failure counts. int runCubicWedgeShapeChecks() { const int kronecker = checkKroneckerDelta(); std::cout << "LinearCubicWedgeShape Kronecker delta failures: " << kronecker diff --git a/test/testLinearWedgeElement.cpp b/test/testLinearWedgeElement.cpp index 7b2638f..0d46d3f 100644 --- a/test/testLinearWedgeElement.cpp +++ b/test/testLinearWedgeElement.cpp @@ -190,4 +190,4 @@ int main(int argc, char **argv) { const int numFailed = runLinearWedgeElementChecks(); Kokkos::finalize(); return numFailed ? 1 : 0; -} +} \ No newline at end of file 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 22/23] 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 From f7daeba18c165ed7f5bfc141cb0fec3702124f6f Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 28 Oct 2026 11:09:21 -0400 Subject: [PATCH 23/23] Pass the controller type to the LinearCubicWedgeShape ShapeField --- src/MeshField_ShapeField.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 7de9091..c7b0b38 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -287,7 +287,7 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { // QuadraticAccessor forwards the dof index to the edge field, so it also // serves edges with more than one dof using QA = QuadraticAccessor; - using WedgeShapeField = ShapeField; + using WedgeShapeField = ShapeField; WedgeShapeField wsf(meshInfo, {vtxField, edgeField}); return FieldWithController{kk_ctrl, wsf}; } else if constexpr (order == 1 && (dim == 1 || dim == 2 || dim == 3)) {