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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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/45] 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 f49d13e2a13244f2e7dae27ca655b06b4f074ad3 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 14:43:51 -0400 Subject: [PATCH 19/45] add a degree-4 triangle integration rule --- src/MeshField_Integrate.hpp | 24 +++++++++++++++++++++++- 1 file changed, 23 insertions(+), 1 deletion(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index da80ba4..56164ef 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -74,11 +74,33 @@ class TriangleIntegration : public EntityIntegration<2> { } virtual int getAccuracy() const { return 2; } }; // end N2 + + // 6-point rule with positive weights, exact for polynomials of degree 4 + // (D. A. Dunavant, Int. J. Numer. Meth. Eng. 21, 1985). The points are the + // two orbits of the barycentric point (a, a, 1-2a). It also serves order 3: + // the only symmetric 4-point degree-3 rule has a negative weight. + class N4 : public Integration<2> { + public: + virtual int countPoints() const { return 6; } + virtual std::vector> getPoints() const { + const Real a1 = 0.44594849091596488632, w1 = 0.22338158967801146570; + const Real a2 = 0.091576213509770743460, w2 = 0.10995174365532186764; + return {IntegrationPoint(Vector2{a1, a1}, w1 / 2.0, 2, 0), + IntegrationPoint(Vector2{1 - 2 * a1, a1}, w1 / 2.0, 2, 0), + IntegrationPoint(Vector2{a1, 1 - 2 * a1}, w1 / 2.0, 2, 0), + IntegrationPoint(Vector2{a2, a2}, w2 / 2.0, 2, 0), + IntegrationPoint(Vector2{1 - 2 * a2, a2}, w2 / 2.0, 2, 0), + IntegrationPoint(Vector2{a2, 1 - 2 * a2}, w2 / 2.0, 2, 0)}; + } + virtual int getAccuracy() const { return 4; } + }; // end N4 + virtual int countIntegrations() const { return 2; } virtual Integration<2> const *getIntegration(int i) const { static N1 i1; static N2 i2; - static Integration<2> *integrations[2] = {&i1, &i2}; + static N4 i4; + static Integration<2> *integrations[3] = {&i1, &i2, &i4}; return integrations[i]; } }; From a1c09552fb92f44e57a2b45caee42f3852968c3d Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 15:00:22 -0400 Subject: [PATCH 20/45] add TriangleIntegrationTest for the triangle quadrature rules --- test/testTriangleIntegration.cpp | 139 +++++++++++++++++++++++++++++++ 1 file changed, 139 insertions(+) create mode 100644 test/testTriangleIntegration.cpp diff --git a/test/testTriangleIntegration.cpp b/test/testTriangleIntegration.cpp new file mode 100644 index 0000000..fd6c6b9 --- /dev/null +++ b/test/testTriangleIntegration.cpp @@ -0,0 +1,139 @@ +#include "KokkosController.hpp" +#include "MeshField.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_Integrate.hpp" +#include "Omega_h_build.hpp" +#include +#include +#include + +// Checks MeshField::TriangleIntegration: +// 1. every rule has positive weights and points inside the reference +// triangle, and integrates xi0^a xi1^b exactly for a + b up to its +// accuracy, +// 2. getAccurate(order) picks the expected rule for orders 1 to 4 and none +// for order 5, +// 3. Integrator::process with orders 3 and 4 integrates monomials of that +// degree exactly over a triangulated unit square. +// Failures are counted explicitly so the test also fails in builds where +// assert() is disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; + +bool isClose(Real a, Real b, Real tol) { + return std::fabs(a - b) <= tol * (1 + std::fabs(b)); +} + +Real factorial(int n) { return n <= 1 ? 1.0 : n * factorial(n - 1); } + +// exact integral of xi0^a xi1^b over the reference triangle +Real triangleMonomial(int a, int b) { + return factorial(a) * factorial(b) / factorial(a + b + 2); +} + +// A rule of accuracy p must integrate xi0^a xi1^b exactly for a + b <= p. +// The a = b = 0 case checks that the weights sum to 1/2, the area of the +// reference triangle. +int checkRule(MeshField::Integration<2> const &rule) { + const auto pts = rule.getPoints(); + const int p = rule.getAccuracy(); + int failed = (static_cast(pts.size()) != rule.countPoints()); + for (auto const &pt : pts) { + const Real x0 = pt.param[0], x1 = pt.param[1]; + failed += !(pt.weight > 0 && x0 >= 0 && x1 >= 0 && x0 + x1 <= 1); + } + for (int a = 0; a <= p; ++a) { + for (int b = 0; a + b <= p; ++b) { + Real sum = 0; + for (auto const &pt : pts) + sum += pt.weight * std::pow(pt.param[0], a) * std::pow(pt.param[1], b); + failed += !isClose(sum, triangleMonomial(a, b), 1e-14); + } + } + std::cout << "triangle rule, accuracy " << p << ", " << pts.size() + << " points: failures " << failed << "\n"; + return failed; +} + +// getAccurate(order) returns the first rule with accuracy >= order +int checkRuleSelection() { + const auto tri = MeshField::getIntegration(); + const int numPts[] = {1, 3, 6, 6}; // orders 1 to 4 + int failed = 0; + for (int order = 1; order <= 4; ++order) + failed += tri->getAccurate(order)->countPoints() != numPts[order - 1]; + failed += tri->getAccurate(5) != nullptr; + std::cout << "rule selection: failures " << failed << "\n"; + return failed; +} + +// Integrates x^px y^py over the mesh. The coordinate map is affine on each +// triangle, so the integrand has degree px + py in the parametric +// coordinates and a rule of that order is exact. +template +class MonomialIntegrator : public MeshField::Integrator { +public: + MonomialIntegrator(FieldElement &feIn, int order, int pxIn, int pyIn) + : MeshField::Integrator(order), fe(feIn), px(pxIn), py(pyIn) {} + void atPoints(Kokkos::View p, Kokkos::View w, + Kokkos::View dV) override { + const int numPtsPerElem = p.extent(0) / fe.numMeshEnts; + const auto x = MeshField::evaluate(fe, p, numPtsPerElem); + const Real a = px, b = py; + Real sum = 0; + Kokkos::parallel_reduce( + "integrateMonomial", w.extent(0), + KOKKOS_LAMBDA(const int i, Real &lsum) { + lsum += + w(i) * dV(i) * Kokkos::pow(x(i, 0), a) * Kokkos::pow(x(i, 1), b); + }, + sum); + integral += sum; + } + FieldElement &fe; + int px, py; + Real integral = 0; +}; + +// integral of x^px y^py over the unit square +Real squareMonomial(int px, int py) { return 1.0 / ((px + 1) * (py + 1)); } + +int checkIntegrator(Omega_h::Mesh &mesh) { + MeshField::OmegahMeshField omf(mesh); + auto coordField = omf.getCoordField(); + const auto [shp, map] = MeshField::Omegah::getTriangleElement<1>(mesh); + MeshField::FieldElement fe(mesh.nelems(), coordField.field, shp, map); + // {order, px, py}: each monomial has degree equal to the order + const int cases[][3] = {{3, 2, 1}, {3, 0, 3}, {4, 2, 2}, {4, 3, 1}}; + int failed = 0; + std::cout.precision(16); + for (auto const &c : cases) { + MonomialIntegrator integrator(fe, c[0], c[1], c[2]); + integrator.process(fe); + const Real expected = squareMonomial(c[1], c[2]); + const bool ok = isClose(integrator.integral, expected, 1e-13); + std::cout << "order " << c[0] << ": integral of x^" << c[1] << " y^" << c[2] + << " = " << integrator.integral << ", expected " << expected + << (ok ? "" : " FAIL") << "\n"; + failed += !ok; + } + return failed; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + int failed = 0; + { + auto lib = Omega_h::Library(&argc, &argv); + failed += checkRule(MeshField::TriangleIntegration::N1()); + failed += checkRule(MeshField::TriangleIntegration::N2()); + failed += checkRule(MeshField::TriangleIntegration::N4()); + failed += checkRuleSelection(); + auto tris = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, + 3, 3, 0); + failed += checkIntegrator(tris); + } + Kokkos::finalize(); + return failed ? 1 : 0; +} From 9af321daafcdcd32dfc964685ce2575b673faaf0 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 15:01:03 -0400 Subject: [PATCH 21/45] remove comment --- test/testTriangleIntegration.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/test/testTriangleIntegration.cpp b/test/testTriangleIntegration.cpp index fd6c6b9..5749068 100644 --- a/test/testTriangleIntegration.cpp +++ b/test/testTriangleIntegration.cpp @@ -15,8 +15,6 @@ // for order 5, // 3. Integrator::process with orders 3 and 4 integrates monomials of that // degree exactly over a triangulated unit square. -// Failures are counted explicitly so the test also fails in builds where -// assert() is disabled. using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; From 8205b1a53a73de60f98a01a073b213f25cce78fc Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 15:20:32 -0400 Subject: [PATCH 22/45] Register TriangleIntegrationTest in CMake --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 45dc805..e8efcf2 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -213,6 +213,7 @@ 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(TriangleIntegrationTest test/testTriangleIntegration.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) if(MeshFields_USE_Cabana) @@ -235,6 +236,7 @@ test_func(CountIntegrator ./CountIntegrator) smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(OmegahTetTest ./OmegahTetTest) test_func(IntegratorMeasure ./IntegratorMeasure) +test_func(TriangleIntegrationTest ./TriangleIntegrationTest) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) From 995db181f36093acb9cccd0b08db5e37b600a413 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 22 Oct 2026 15:40:49 -0400 Subject: [PATCH 23/45] set countIntegrations() to return 3 --- src/MeshField_Integrate.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index 56164ef..5dcb2b1 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -95,7 +95,7 @@ class TriangleIntegration : public EntityIntegration<2> { virtual int getAccuracy() const { return 4; } }; // end N4 - virtual int countIntegrations() const { return 2; } + virtual int countIntegrations() const { return 3; } virtual Integration<2> const *getIntegration(int i) const { static N1 i1; static N2 i2; From 59e41b120035cd338c58d778a56e823380b64556 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Mon, 26 Oct 2026 14:30:42 -0400 Subject: [PATCH 24/45] add reference --- src/MeshField_Integrate.hpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index 5dcb2b1..f425889 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -75,10 +75,13 @@ class TriangleIntegration : public EntityIntegration<2> { virtual int getAccuracy() const { return 2; } }; // end N2 - // 6-point rule with positive weights, exact for polynomials of degree 4 - // (D. A. Dunavant, Int. J. Numer. Meth. Eng. 21, 1985). The points are the + // 6-point rule with positive weights, exact for polynomials of degree 4. The points are the // two orbits of the barycentric point (a, a, 1-2a). It also serves order 3: // the only symmetric 4-point degree-3 rule has a negative weight. + // + // D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature + // rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). + // class N4 : public Integration<2> { public: virtual int countPoints() const { return 6; } From e542ba04cbc5d2729386199186e1b6c7afe53666 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Mon, 26 Oct 2026 15:26:30 -0400 Subject: [PATCH 25/45] copied N4 rule from apfIntegrate.cc --- src/MeshField_Integrate.hpp | 36 ++++++++++++++++++++---------------- 1 file changed, 20 insertions(+), 16 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index f425889..c2abdec 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -75,29 +75,33 @@ class TriangleIntegration : public EntityIntegration<2> { virtual int getAccuracy() const { return 2; } }; // end N2 - // 6-point rule with positive weights, exact for polynomials of degree 4. The points are the - // two orbits of the barycentric point (a, a, 1-2a). It also serves order 3: - // the only symmetric 4-point degree-3 rule has a negative weight. - // - // D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature - // rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). - // +// Symmetric Gaussian quadrature rules for the triangle from +// D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature +// rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). +// N4 is p = 4 rules of from the Appendix II (p. 1140), +// copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). +// weights sum to one, so divided by 2, the area of the reference +// triangle. The p = 3 rule has a negative weight, so order 3 uses N4. class N4 : public Integration<2> { public: virtual int countPoints() const { return 6; } virtual std::vector> getPoints() const { - const Real a1 = 0.44594849091596488632, w1 = 0.22338158967801146570; - const Real a2 = 0.091576213509770743460, w2 = 0.10995174365532186764; - return {IntegrationPoint(Vector2{a1, a1}, w1 / 2.0, 2, 0), - IntegrationPoint(Vector2{1 - 2 * a1, a1}, w1 / 2.0, 2, 0), - IntegrationPoint(Vector2{a1, 1 - 2 * a1}, w1 / 2.0, 2, 0), - IntegrationPoint(Vector2{a2, a2}, w2 / 2.0, 2, 0), - IntegrationPoint(Vector2{1 - 2 * a2, a2}, w2 / 2.0, 2, 0), - IntegrationPoint(Vector2{a2, 1 - 2 * a2}, w2 / 2.0, 2, 0)}; + return {IntegrationPoint(Vector2{0.816847572980459, 0.091576213509771}, + 0.109951743655322 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.091576213509771, 0.816847572980459}, + 0.109951743655322 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.091576213509771, 0.091576213509771}, + 0.109951743655322 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.108103018168070, 0.445948490915965}, + 0.223381589678011 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.445948490915965, 0.108103018168070}, + 0.223381589678011 / 2.0, 2, 0), + IntegrationPoint(Vector2{0.445948490915965, 0.445948490915965}, + 0.223381589678011 / 2.0, 2, 0)}; } virtual int getAccuracy() const { return 4; } }; // end N4 - + virtual int countIntegrations() const { return 3; } virtual Integration<2> const *getIntegration(int i) const { static N1 i1; 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 26/45] 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 27/45] 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 93cf6d7bc0af8de5e7e92ba5822cd61e37179d74 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 00:41:00 -0400 Subject: [PATCH 28/45] add wedge integration points and weights --- src/MeshField_Integrate.hpp | 177 ++++++++++++++++++++++++++++++++---- 1 file changed, 159 insertions(+), 18 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index c2abdec..c1a8a15 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -47,6 +47,101 @@ template class EntityIntegration { virtual int countIntegrations() const = 0; virtual Integration const *getIntegration(int i) const = 0; }; +// Gauss-Legendre rules on [-1, 1] from +// N1 to N6 are copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc +// (EdgeIntegration, lines 27-129). The n-point rule is exact for degree +// 2n - 1, and the weights sum to 2, the length of the reference edge. +//! [EdgeIntegration] +class EdgeIntegration : public EntityIntegration<1> { +public: + class N1 : public Integration<1> { + public: + virtual int countPoints() const { return 1; } + virtual std::vector> getPoints() const { + return {IntegrationPoint(Vector1{0}, 2, 1, 0)}; + } + virtual int getAccuracy() const { return 1; } + }; // end N1 + class N2 : public Integration<1> { + public: + virtual int countPoints() const { return 2; } + virtual std::vector> getPoints() const { + return {IntegrationPoint(Vector1{0.577350269189626}, 1, 1, 0), + IntegrationPoint(Vector1{-0.577350269189626}, 1, 1, 0)}; + } + virtual int getAccuracy() const { return 3; } + }; // end N2 + class N3 : public Integration<1> { + public: + virtual int countPoints() const { return 3; } + virtual std::vector> getPoints() const { + return { + IntegrationPoint(Vector1{0.000000000000000}, 0.888888888888889, 1, 0), + IntegrationPoint(Vector1{0.774596669241483}, 0.555555555555556, 1, 0), + IntegrationPoint(Vector1{-0.774596669241483}, 0.555555555555556, 1, + 0)}; + } + virtual int getAccuracy() const { return 5; } + }; // end N3 + class N4 : public Integration<1> { + public: + virtual int countPoints() const { return 4; } + virtual std::vector> getPoints() const { + return { + IntegrationPoint(Vector1{0.339981043584856}, 0.652145154862546, 1, 0), + IntegrationPoint(Vector1{-0.339981043584856}, 0.652145154862546, 1, + 0), + IntegrationPoint(Vector1{0.861136311594053}, 0.347854845137454, 1, 0), + IntegrationPoint(Vector1{-0.861136311594053}, 0.347854845137454, 1, + 0)}; + } + virtual int getAccuracy() const { return 7; } + }; // end N4 + class N5 : public Integration<1> { + public: + virtual int countPoints() const { return 5; } + virtual std::vector> getPoints() const { + return { + IntegrationPoint(Vector1{0.000000000000000}, 0.568888888888889, 1, 0), + IntegrationPoint(Vector1{0.538469310105683}, 0.478628670499366, 1, 0), + IntegrationPoint(Vector1{-0.538469310105683}, 0.478628670499366, 1, + 0), + IntegrationPoint(Vector1{0.906179845938664}, 0.236926885056189, 1, 0), + IntegrationPoint(Vector1{-0.906179845938664}, 0.236926885056189, 1, + 0)}; + } + virtual int getAccuracy() const { return 9; } + }; // end N5 + class N6 : public Integration<1> { + public: + virtual int countPoints() const { return 6; } + virtual std::vector> getPoints() const { + return { + IntegrationPoint(Vector1{0.238619186083197}, 0.467913934572691, 1, 0), + IntegrationPoint(Vector1{-0.238619186083197}, 0.467913934572691, 1, + 0), + IntegrationPoint(Vector1{0.661209386466265}, 0.360761573048139, 1, 0), + IntegrationPoint(Vector1{-0.661209386466265}, 0.360761573048139, 1, + 0), + IntegrationPoint(Vector1{0.932469514203152}, 0.171324492379170, 1, 0), + IntegrationPoint(Vector1{-0.932469514203152}, 0.171324492379170, 1, + 0)}; + } + virtual int getAccuracy() const { return 11; } + }; // end N6 + virtual int countIntegrations() const { return 6; } + virtual Integration<1> const *getIntegration(int i) const { + static N1 i1; + static N2 i2; + static N3 i3; + static N4 i4; + static N5 i5; + static N6 i6; + static Integration<1> *integrations[6] = {&i1, &i2, &i3, &i4, &i5, &i6}; + return integrations[i]; + } +}; +//! [EdgeIntegration] // points are given in the reduced parametric coordinates (xi0=L1, xi1=L2, // ...); the redundant barycentric coordinate L0 = 1-sum(xi) is omitted @@ -65,23 +160,23 @@ class TriangleIntegration : public EntityIntegration<2> { public: virtual int countPoints() const { return 3; } virtual std::vector> getPoints() const { - return {IntegrationPoint( - Vector2{0.166666666666667, 0.166666666666667}, 1. / 3. / 2.0, 2, 0), - IntegrationPoint( - Vector2{0.666666666666667, 0.166666666666667}, 1. / 3. / 2.0, 2, 0), - IntegrationPoint( - Vector2{0.166666666666667, 0.666666666666667}, 1. / 3. / 2.0, 2, 0)}; + return {IntegrationPoint(Vector2{0.166666666666667, 0.166666666666667}, + 1. / 3. / 2.0, 2, 0), + IntegrationPoint(Vector2{0.666666666666667, 0.166666666666667}, + 1. / 3. / 2.0, 2, 0), + IntegrationPoint(Vector2{0.166666666666667, 0.666666666666667}, + 1. / 3. / 2.0, 2, 0)}; } virtual int getAccuracy() const { return 2; } }; // end N2 -// Symmetric Gaussian quadrature rules for the triangle from -// D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature -// rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). -// N4 is p = 4 rules of from the Appendix II (p. 1140), -// copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). -// weights sum to one, so divided by 2, the area of the reference -// triangle. The p = 3 rule has a negative weight, so order 3 uses N4. + // Symmetric Gaussian quadrature rules for the triangle from + // D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature + // rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). + // N4 is p = 4 rules of from the Appendix II (p. 1140), + // copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). + // weights sum to one, so divided by 2, the area of the reference + // triangle. The p = 3 rule has a negative weight, so order 3 uses N4. class N4 : public Integration<2> { public: virtual int countPoints() const { return 6; } @@ -101,7 +196,7 @@ class TriangleIntegration : public EntityIntegration<2> { } virtual int getAccuracy() const { return 4; } }; // end N4 - + virtual int countIntegrations() const { return 3; } virtual Integration<2> const *getIntegration(int i) const { static N1 i1; @@ -113,6 +208,7 @@ class TriangleIntegration : public EntityIntegration<2> { }; //! [TriangleIntegration] +//! [TetrahedronIntegrationIntegration] class TetrahedronIntegration : public EntityIntegration<3> { public: class N1 : public Integration<3> { @@ -151,19 +247,52 @@ class TetrahedronIntegration : public EntityIntegration<3> { return integrations[i]; } }; +//! [TetrahedronIntegrationIntegration] + +// Wedge rules are the tensor product of a triangle rule in (xi0, xi1) and an +// edge rule in zeta. The rule for (triOrder, zetaOrder) integrates +// xi0^a xi1^b zeta^c exactly for a + b <= triOrder and c <= zetaOrder. The +// weights sum to 1, the volume of the reference wedge (area 1/2 times length +// 2). +inline std::vector> +getWedgeIntegrationPoints(int triOrder, int zetaOrder) { + const auto *tri = TriangleIntegration().getAccurate(triOrder); + if (tri == nullptr) { + fail("getWedgeIntegrationPoints: no triangle rule of order %d\n", triOrder); + } + const auto *edge = EdgeIntegration().getAccurate(zetaOrder); + if (edge == nullptr) { + fail("getWedgeIntegrationPoints: no edge rule of order %d\n", zetaOrder); + } + std::vector> points; + for (auto const &t : tri->getPoints()) + for (auto const &z : edge->getPoints()) + points.push_back( + IntegrationPoint(Vector3{t.param[0], t.param[1], z.param[0]}, + t.weight * z.weight, 3, 0)); + return points; +} + //! [getIntegration] template auto const getIntegration() { - if constexpr (topo == Triangle) { + if constexpr (topo == Edge) { + return std::make_shared(); + } else if constexpr (topo == Triangle) { return std::make_shared(); } else if constexpr (topo == Tetrahedron) { return std::make_shared(); } fail("getIntegration does not support given topology\n"); } + //! [getIntegration] template auto getIntegrationPoints(int order) { - auto ip = getIntegration()->getAccurate(order)->getPoints(); - return ip; + if constexpr (topo == Wedge) { + return getWedgeIntegrationPoints(order, order); + } else { + auto ip = getIntegration()->getAccurate(order)->getPoints(); + return ip; + } } template @@ -245,6 +374,11 @@ class Integrator { public: /** \brief Construct an Integrator given an order of accuracy. */ Integrator(int o) : order(o) {} + /** \brief Construct an Integrator for wedges, with separate orders of + * accuracy for the triangle (xi0, xi1) and zeta directions. */ + Integrator(int triOrder, int zetaOrderIn) + : order(triOrder), zetaOrder(zetaOrderIn) {} + virtual ~Integrator(){}; /** \brief User callback: process entry * @@ -284,7 +418,13 @@ class Integrator { template void process(FieldElement &fes) { constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology(); pre(); - auto ip = getIntegrationPoints(order); + auto ip = [&] { + if constexpr (topo[0] == Wedge) { + return getWedgeIntegrationPoints(order, zetaOrder); + } else { + return getIntegrationPoints(order); + } + }(); auto localCoords = getIntegrationPointLocalCoords(fes, ip); auto weights = getIntegrationPointWeights(fes, ip); auto dV = getJacobianDeterminants(fes, localCoords, ip.size()); @@ -299,6 +439,7 @@ class Integrator { protected: int order; + int zetaOrder; }; } // namespace MeshField From 65b1ce78fd88058c2d177f632bbb9eaa6ec91ad9 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 00:47:24 -0400 Subject: [PATCH 29/45] add test case to test edge quadrature rule and wedge tensor prdouct rule --- test/testWedgeIntegration.cpp | 196 ++++++++++++++++++++++++++++++++++ 1 file changed, 196 insertions(+) create mode 100644 test/testWedgeIntegration.cpp diff --git a/test/testWedgeIntegration.cpp b/test/testWedgeIntegration.cpp new file mode 100644 index 0000000..65b36d2 --- /dev/null +++ b/test/testWedgeIntegration.cpp @@ -0,0 +1,196 @@ +#include "KokkosController.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_For.hpp" +#include "MeshField_Integrate.hpp" +#include "MeshField_ShapeField.hpp" +#include +#include +#include + +// Checks the wedge quadrature: +// 1. every EdgeIntegration rule has positive weights and points inside +// (-1, 1), and integrates zeta^c exactly for c up to its accuracy, +// 2. getWedgeIntegrationPoints(triOrder, zetaOrder) has the expected number +// of points, positive weights and points inside the reference wedge, and +// integrates xi0^a xi1^b zeta^c exactly for a + b <= triOrder and +// c <= zetaOrder, +// 3. Integrator::process with LinearWedgeShape gives the exact volume of a +// right prism and a tapered wedge, and a zeta order that is too low gives +// the wrong volume for the tapered wedge. +// Failures are counted explicitly so the test also fails in builds where +// assert() is disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; + +bool isClose(Real a, Real b, Real tol) { + return std::fabs(a - b) <= tol * (1 + std::fabs(b)); +} + +Real factorial(int n) { return n <= 1 ? 1.0 : n * factorial(n - 1); } + +// exact integral of xi0^a xi1^b over the reference triangle +Real triangleMonomial(int a, int b) { + return factorial(a) * factorial(b) / factorial(a + b + 2); +} + +// exact integral of zeta^c over [-1, 1] +Real lineMonomial(int c) { return (c % 2 == 1) ? 0.0 : 2.0 / (c + 1); } + +// A rule of accuracy p must integrate zeta^c exactly for c <= p. The c = 0 +// case checks that the weights sum to 2, the length of the reference edge. +int checkEdgeRule(MeshField::Integration<1> const &rule) { + const auto pts = rule.getPoints(); + const int p = rule.getAccuracy(); + int failed = (static_cast(pts.size()) != rule.countPoints()); + for (auto const &pt : pts) + failed += !(pt.weight > 0 && pt.param[0] > -1 && pt.param[0] < 1); + for (int c = 0; c <= p; ++c) { + Real sum = 0; + for (auto const &pt : pts) + sum += pt.weight * std::pow(pt.param[0], c); + failed += !isClose(sum, lineMonomial(c), 1e-14); + } + std::cout << "edge rule, accuracy " << p << ", " << pts.size() + << " points: failures " << failed << "\n"; + return failed; +} + +int checkWedgeRule(int triOrder, int zetaOrder, int numTriPts, + int numGaussPts) { + const auto pts = MeshField::getWedgeIntegrationPoints(triOrder, zetaOrder); + int failed = (static_cast(pts.size()) != numTriPts * numGaussPts); + for (auto const &pt : pts) { + const Real x0 = pt.param[0], x1 = pt.param[1], z = pt.param[2]; + failed += !(pt.weight > 0 && x0 >= 0 && x1 >= 0 && x0 + x1 <= 1 && z > -1 && + z < 1); + } + // a = b = c = 0 checks that the weights sum to 1 + for (int a = 0; a <= triOrder; ++a) + for (int b = 0; a + b <= triOrder; ++b) + for (int c = 0; c <= zetaOrder; ++c) { + Real sum = 0; + for (auto const &pt : pts) + sum += pt.weight * std::pow(pt.param[0], a) * + std::pow(pt.param[1], b) * std::pow(pt.param[2], c); + failed += + !isClose(sum, triangleMonomial(a, b) * lineMonomial(c), 1e-14); + } + std::cout << "wedge rule (" << triOrder << ", " << zetaOrder << "), " + << pts.size() << " points: failures " << failed << "\n"; + return failed; +} + +// node n of wedge w is vertex 6*w + n +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); + return {0, wedgeCompIdx, 6 * wedge + wedgeNodeIdx, MeshField::Vertex}; + } +}; + +// sums w * |det J| over every point of every wedge +class VolumeIntegrator : public MeshField::Integrator { +public: + VolumeIntegrator(int triOrder, int zetaOrder) + : MeshField::Integrator(triOrder, zetaOrder) {} + void atPoints(Kokkos::View, Kokkos::View w, + Kokkos::View dV) override { + Real sum = 0; + Kokkos::parallel_reduce( + "volume", w.extent(0), + KOKKOS_LAMBDA(const int i, Real &s) { s += w(i) * dV(i); }, sum); + volume = sum; + } + Real volume = 0; +}; + +// clang-format off +// wedge 0: right prism, triangle (0,0),(2,0),(0,3) times height 4, volume 12 +// wedge 1: tapered wedge, bottom triangle (0,0),(2,0),(0,2) at z = 0 and top +// triangle (0,0),(1,0),(0,1) at z = 1, a frustum with volume +// h/3 (A0 + A1 + sqrt(A0 A1)) = (2 + 1/2 + 1)/3 = 7/6 +constexpr int numWedges = 2; +const Real wedgeCoords[numWedges][6][3] = { + {{0, 0, 0}, {2, 0, 0}, {0, 3, 0}, {0, 0, 4}, {2, 0, 4}, {0, 3, 4}}, + {{0, 0, 0}, {2, 0, 0}, {0, 2, 0}, {0, 0, 1}, {1, 0, 1}, {0, 1, 1}}}; +// clang-format on +const Real exactVolume = 12.0 + 7.0 / 6.0; + +template void setVtxCoords(CoordField field) { + const int numVtx = 6 * numWedges; + Kokkos::View coords_h("coords_h", numVtx); + for (int w = 0; w < numWedges; ++w) + for (int n = 0; n < 6; ++n) + for (int d = 0; d < 3; ++d) + coords_h(6 * 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"); +} + +int checkIntegrator() { + const MeshField::MeshInfo meshInfo{ + .numVtx = 6 * numWedges, .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()); + // det J of a linear wedge is degree 1 in (xi0, xi1) and degree 2 in zeta + const int cases[][2] = {{1, 2}, {2, 2}, {3, 11}}; + int failed = 0; + std::cout.precision(16); + for (auto const &c : cases) { + VolumeIntegrator vi(c[0], c[1]); + vi.process(fe); + const bool ok = isClose(vi.volume, exactVolume, 1e-13); + failed += !ok; + std::cout << "Integrator (" << c[0] << ", " << c[1] << "): volume " + << vi.volume << ", expected " << exactVolume + << (ok ? "" : " FAIL") << "\n"; + } + // one Gauss point in zeta is not enough for the tapered wedge + VolumeIntegrator low(1, 1); + low.process(fe); + const bool lowDiffers = !isClose(low.volume, exactVolume, 1e-6); + failed += !lowDiffers; + std::cout << "Integrator (1, 1): volume " << low.volume + << (lowDiffers ? " (not exact, as expected)" : " FAIL") << "\n"; + return failed; +} + +int runWedgeIntegrationChecks() { + int failed = checkEdgeRule(MeshField::EdgeIntegration::N1()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N2()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N3()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N4()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N5()); + failed += checkEdgeRule(MeshField::EdgeIntegration::N6()); + failed += checkWedgeRule(1, 1, 1, 1); + failed += checkWedgeRule(2, 2, 3, 2); + failed += checkWedgeRule(2, 5, 3, 3); + failed += checkWedgeRule(4, 4, 6, 3); + failed += checkWedgeRule(3, 11, 6, 6); + failed += checkIntegrator(); + return failed; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + const int numFailed = runWedgeIntegrationChecks(); + Kokkos::finalize(); + return numFailed ? 1 : 0; +} From 4df32bc5728500c4d40fcace92c83c8471bf90ff Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 00:49:43 -0400 Subject: [PATCH 30/45] register test case in cmakelists --- CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/CMakeLists.txt b/CMakeLists.txt index 415ef1d..d9ec59f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -216,6 +216,7 @@ meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp) meshfields_add_exe(LinearWedgeShapeTest test/testLinearWedgeShape.cpp) meshfields_add_exe(LinearWedgeElementTest test/testLinearWedgeElement.cpp) meshfields_add_exe(TriangleIntegrationTest test/testTriangleIntegration.cpp) +meshfields_add_exe(WedgeIntegrationTest test/testWedgeIntegration.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) if(MeshFields_USE_Cabana) @@ -241,6 +242,7 @@ test_func(IntegratorMeasure ./IntegratorMeasure) test_func(LinearWedgeShapeTest ./LinearWedgeShapeTest) test_func(LinearWedgeElementTest ./LinearWedgeElementTest) test_func(TriangleIntegrationTest ./TriangleIntegrationTest) +test_func(WedgeIntegrationTest ./WedgeIntegrationTest) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) From ff385214e7dd6b902e7211efeef781ddc36115fd Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 01:07:26 -0400 Subject: [PATCH 31/45] call getIntegrationPoints directly in Integrator::process --- src/MeshField_Integrate.hpp | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index c1a8a15..43fb0fd 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -284,14 +284,16 @@ template auto const getIntegration() { } fail("getIntegration does not support given topology\n"); } - //! [getIntegration] -template auto getIntegrationPoints(int order) { + +// wedges use triOrder in (xi0, xi1) and zetaOrder in zeta; other topologies +// use triOrder +template +auto getIntegrationPoints(int triOrder, int zetaOrder) { if constexpr (topo == Wedge) { - return getWedgeIntegrationPoints(order, order); + return getWedgeIntegrationPoints(triOrder, zetaOrder); } else { - auto ip = getIntegration()->getAccurate(order)->getPoints(); - return ip; + return getIntegrationPoints(triOrder); } } @@ -418,13 +420,7 @@ class Integrator { template void process(FieldElement &fes) { constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology(); pre(); - auto ip = [&] { - if constexpr (topo[0] == Wedge) { - return getWedgeIntegrationPoints(order, zetaOrder); - } else { - return getIntegrationPoints(order); - } - }(); + auto ip = getIntegrationPoints(order, zetaOrder); auto localCoords = getIntegrationPointLocalCoords(fes, ip); auto weights = getIntegrationPointWeights(fes, ip); auto dV = getJacobianDeterminants(fes, localCoords, ip.size()); From 33376df7dcf24a90f7106f41c68e1c421f6a9436 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 01:09:25 -0400 Subject: [PATCH 32/45] clang-format the src files --- src/CabanaController.hpp | 4 +- src/MeshField.hpp | 4 +- src/MeshField_Config.hpp.in | 2 +- src/MeshField_Element.hpp | 10 +- src/MeshField_Shape.hpp | 213 ++++++++++++++++++++--------------- src/MeshField_ShapeField.hpp | 14 ++- 6 files changed, 141 insertions(+), 106 deletions(-) diff --git a/src/CabanaController.hpp b/src/CabanaController.hpp index aadf534..289ce53 100644 --- a/src/CabanaController.hpp +++ b/src/CabanaController.hpp @@ -111,8 +111,8 @@ class CabanaController { template using wrapper_slice_t = CabanaSliceWrapper, T>; - using AoSoATuple = - std::tuple, MemorySpace, vecLen>...>; + using AoSoATuple = std::tuple< + Cabana::AoSoA, MemorySpace, vecLen>...>; template static auto construct(std::integer_sequence, diff --git a/src/MeshField.hpp b/src/MeshField.hpp index dd541b8..a4bd8f6 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -223,8 +223,8 @@ struct QuadraticTetrahedronToField { const auto tetDim = 3; const auto vtxDim = 0; const auto ignored = -1; - // cyclic rotation of the omegah vertex order to map to the meshfields order - // defined by the shape functions in MeshField_Shape.hpp + // cyclic rotation of the omegah vertex order to map to the meshfields + // order defined by the shape functions in MeshField_Shape.hpp const auto localVtxIdx = (Omega_h::simplex_down_template( tetDim, vtxDim, dofHolderIdx, ignored) + 3) % diff --git a/src/MeshField_Config.hpp.in b/src/MeshField_Config.hpp.in index cb05f97..5287cd1 100644 --- a/src/MeshField_Config.hpp.in +++ b/src/MeshField_Config.hpp.in @@ -1,2 +1,2 @@ #cmakedefine MeshFields_USE_EXCEPTIONS -#cmakedefine MeshFields_VERSION @MeshFields_VERSION@ +#cmakedefine MeshFields_VERSION @MeshFields_VERSION @ diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index cb2890d..6ed5b1b 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -88,17 +88,19 @@ KOKKOS_INLINE_FUNCTION auto addTensorProduct(VecA const &a, VecB const &b, namespace MeshField { - // /** * @brief Supports mapping between mesh (i.e., Omega_h) ordering and MeshFields * ordering * - * Return type used by structs/classes that implement the ElementDofHolderAccessor parenthesis operator. See @ref adding-shape-functions-omegah "Adding Shape Functions Omega_h Element Topologies" + * Return type used by structs/classes that implement the + * ElementDofHolderAccessor parenthesis operator. See @ref + * adding-shape-functions-omegah "Adding Shape Functions Omega_h Element + * Topologies" */ struct ElementToDofHolderMap { - /** Which node is being accessed; non-zero + /** Which node is being accessed; non-zero * for mesh entities associated with multiple nodes */ LO node; @@ -114,7 +116,7 @@ struct ElementToDofHolderMap { */ LO entity; - /** Topological type of the entity at entityIndex + /** Topological type of the entity at entityIndex * (e.g., Vertex, Edge, Triangle) */ Mesh_Topology topo; diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index 8489814..0b0195b 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -23,8 +23,11 @@ namespace { * @param val Reference value * @return true if xi >= val within machine precision */ -KOKKOS_INLINE_FUNCTION bool greaterThanOrEqual(MeshField::Real xi, const MeshField::Real val, const MeshField::Real tol) { - if ( xi > val ) return true; +KOKKOS_INLINE_FUNCTION bool greaterThanOrEqual(MeshField::Real xi, + const MeshField::Real val, + const MeshField::Real tol) { + if (xi > val) + return true; return (Kokkos::fabs(xi - val) <= tol); } @@ -36,10 +39,12 @@ KOKKOS_INLINE_FUNCTION bool greaterThanOrEqual(MeshField::Real xi, const MeshFie * @return true if all xi[i] >= val within machine precision */ template -KOKKOS_INLINE_FUNCTION bool eachGreaterThanOrEqual(Array &xi, const MeshField::Real val, const MeshField::Real tol) { +KOKKOS_INLINE_FUNCTION bool eachGreaterThanOrEqual(Array &xi, + const MeshField::Real val, + const MeshField::Real tol) { auto gt = true; for (size_t i = 0; i < xi.size(); i++) { - gt = gt && greaterThanOrEqual(xi[i],val,tol); + gt = gt && greaterThanOrEqual(xi[i], val, tol); } return gt; } @@ -50,8 +55,11 @@ KOKKOS_INLINE_FUNCTION bool eachGreaterThanOrEqual(Array &xi, const MeshField::R * @param val Reference value * @return true if xi <= val within machine precision */ -KOKKOS_INLINE_FUNCTION bool lessThanOrEqual(MeshField::Real xi, const MeshField::Real val, const MeshField::Real tol) { - if ( xi < val ) return true; +KOKKOS_INLINE_FUNCTION bool lessThanOrEqual(MeshField::Real xi, + const MeshField::Real val, + const MeshField::Real tol) { + if (xi < val) + return true; return (Kokkos::fabs(xi - val) <= tol); } @@ -63,10 +71,12 @@ KOKKOS_INLINE_FUNCTION bool lessThanOrEqual(MeshField::Real xi, const MeshField: * @return true if all xi[i] <= val within machine precision */ template -KOKKOS_INLINE_FUNCTION bool eachLessThanOrEqual(Array &xi, const MeshField::Real val, const MeshField::Real tol) { +KOKKOS_INLINE_FUNCTION bool eachLessThanOrEqual(Array &xi, + const MeshField::Real val, + const MeshField::Real tol) { auto lt = true; for (size_t i = 0; i < xi.size(); i++) { - lt = lt && lessThanOrEqual(xi[i],val,tol); + lt = lt && lessThanOrEqual(xi[i], val, tol); } return lt; } @@ -97,10 +107,11 @@ using Vector3 = Kokkos::Array; * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct LinearEdgeShape { - static const size_t numNodes = 2; ///< Number of nodes (2 for linear edge) - static const size_t meshEntDim = 1; ///< Mesh entity dimension (1D) - constexpr static Mesh_Topology DofHolders[1] = {Vertex}; ///< DOFs located at vertices - constexpr static size_t Order = 1; ///< Polynomial order (linear) + static const size_t numNodes = 2; ///< Number of nodes (2 for linear edge) + static const size_t meshEntDim = 1; ///< Mesh entity dimension (1D) + constexpr static Mesh_Topology DofHolders[1] = { + Vertex}; ///< DOFs located at vertices + constexpr static size_t Order = 1; ///< Polynomial order (linear) /** * @brief Get parametric coordinates of element nodes @@ -121,8 +132,8 @@ struct LinearEdgeShape { */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector1 const &xi) const { - assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi,-1.0,ParametricCoordTol)); + assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi, -1.0, ParametricCoordTol)); // clang-format off return {(1.0 - xi[0]) / 2.0, (1.0 + xi[0]) / 2.0}; @@ -154,22 +165,26 @@ struct LinearEdgeShape { * - Node 1 at (1, 0) * - Node 2 at (0, 1) * - * Shape functions use barycentric coordinates where L0 = 1-xi0-xi1, L1 = xi0, L2 = xi1 + * Shape functions use barycentric coordinates where L0 = 1-xi0-xi1, L1 = xi0, + * L2 = xi1 * * Shape functions and ordering from: * Zienkiewicz, Taylor, and Zhu * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct LinearTriangleShape { - static const size_t numNodes = 3; ///< Number of nodes (3 for linear triangle) - static const size_t numComponentsPerDof = 1; ///< Components per DOF (scalar field) - static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) - constexpr static Mesh_Topology DofHolders[1] = {Vertex}; ///< DOFs located at vertices - constexpr static size_t Order = 1; ///< Polynomial order (linear) + static const size_t numNodes = 3; ///< Number of nodes (3 for linear triangle) + static const size_t numComponentsPerDof = + 1; ///< Components per DOF (scalar field) + static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) + constexpr static Mesh_Topology DofHolders[1] = { + Vertex}; ///< DOFs located at vertices + constexpr static size_t Order = 1; ///< Polynomial order (linear) /** * @brief Get parametric coordinates of element nodes - * @return Array of node coordinates: [node0_xi0, node0_xi1, node1_xi0, node1_xi1, ...] + * @return Array of node coordinates: [node0_xi0, node0_xi1, node1_xi0, + * node1_xi1, ...] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getNodeParametricCoords() const { @@ -188,11 +203,11 @@ struct LinearTriangleShape { */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector2 const &xi) const { - assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); + assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1]; - assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); - assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); // clang-format off return {L0, xi[0], //L1 @@ -203,11 +218,12 @@ struct LinearTriangleShape { /** * @brief Get gradients of shape functions * @detail gradients are constant - ignoring parametric coord input - * @return Array of gradients [dN0/dxi0, dN0/dxi1, dN1/dxi0, dN1/dxi1, dN2/dxi0, dN2/dxi1] - * Constant gradients: [[-1,-1], [1,0], [0,1]] + * @return Array of gradients [dN0/dxi0, dN0/dxi1, dN1/dxi0, dN1/dxi1, + * dN2/dxi0, dN2/dxi1] Constant gradients: [[-1,-1], [1,0], [0,1]] */ KOKKOS_INLINE_FUNCTION - Kokkos::Array getLocalGradients(Vector2 const &) const { + Kokkos::Array + getLocalGradients(Vector2 const &) const { // clang-format off return { -1,-1, //first vector 1, 0, @@ -228,14 +244,16 @@ struct LinearTriangleShape { * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct LinearTriangleCoordinateShape { - static const size_t numNodes = 3; ///< Number of nodes - static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) - constexpr static Mesh_Topology DofHolders[1] = {Vertex}; ///< DOFs located at vertices - constexpr static size_t Order = 1; ///< Polynomial order (linear) + static const size_t numNodes = 3; ///< Number of nodes + static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) + constexpr static Mesh_Topology DofHolders[1] = { + Vertex}; ///< DOFs located at vertices + constexpr static size_t Order = 1; ///< Polynomial order (linear) /** * @brief Get parametric coordinates of element nodes - * @return Array of node coordinates: [node0_xi0, node0_xi1, node1_xi0, node1_xi1, node2_xi0, node2_xi1] + * @return Array of node coordinates: [node0_xi0, node0_xi1, node1_xi0, + * node1_xi1, node2_xi0, node2_xi1] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getNodeParametricCoords() const { @@ -249,15 +267,16 @@ struct LinearTriangleCoordinateShape { /** * @brief Evaluate shape functions at parametric coordinate * @param xi Parametric coordinates (xi0, xi1) where xi0, xi1 \f$\in\f$ [0,1] - * @return Array of shape function values using barycentric coordinates [L0, L1, L2] + * @return Array of shape function values using barycentric coordinates [L0, + * L1, L2] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector2 const &xi) const { - assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); + assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1]; - assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); - assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); // clang-format off return {L0, xi[0], @@ -283,16 +302,19 @@ struct LinearTriangleCoordinateShape { */ //! [QuadraticTriangleShape] struct QuadraticTriangleShape { - static const size_t numNodes = 6; ///< Number of nodes (6 for quadratic triangle) - static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) - constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge}; ///< DOFs at vertices and edges - constexpr static size_t NumDofHolders[2] = {3, 3}; ///< 3 vertices, 3 edges - constexpr static size_t DofsPerHolder[2] = {1, 1}; ///< 1 DOF per vertex/edge - constexpr static size_t Order = 2; ///< Polynomial order (quadratic) + static const size_t numNodes = + 6; ///< Number of nodes (6 for quadratic triangle) + static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) + constexpr static Mesh_Topology DofHolders[2] = { + Vertex, Edge}; ///< DOFs at vertices and edges + constexpr static size_t NumDofHolders[2] = {3, 3}; ///< 3 vertices, 3 edges + constexpr static size_t DofsPerHolder[2] = {1, 1}; ///< 1 DOF per vertex/edge + constexpr static size_t Order = 2; ///< Polynomial order (quadratic) /** * @brief Get parametric coordinates of element nodes - * @return Array of node coordinates for 6 nodes (3 at vertices, 3 at edge midpoints) + * @return Array of node coordinates for 6 nodes (3 at vertices, 3 at edge + * midpoints) */ KOKKOS_INLINE_FUNCTION Kokkos::Array getNodeParametricCoords() const { @@ -320,11 +342,11 @@ struct QuadraticTriangleShape { */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector2 const &xi) const { - assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); + assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1]; - assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); - assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); const Real L1 = xi[0]; const Real L2 = xi[1]; // clang-format off @@ -340,16 +362,17 @@ struct QuadraticTriangleShape { /** * @brief Get gradients of quadratic shape functions at parametric coordinate * @param xi Parametric coordinates (xi0, xi1) where xi0, xi1 \f$\in\f$ [0,1] - * @return Array of gradient vectors [dN0/dxi0, dN0/dxi1, ..., dN5/dxi0, dN5/dxi1] + * @return Array of gradient vectors [dN0/dxi0, dN0/dxi1, ..., dN5/dxi0, + * dN5/dxi1] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getLocalGradients(Vector2 const &xi) const { - assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); + assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1]; - assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); - assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); const Real L1 = xi[0]; const Real L2 = xi[1]; // clang-format off @@ -386,10 +409,12 @@ struct QuadraticTriangleShape { * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct LinearTetrahedronShape { - static const size_t numNodes = 4; ///< Number of nodes (4 for linear tetrahedron) - 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) + static const size_t numNodes = + 4; ///< Number of nodes (4 for linear tetrahedron) + 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) /** * @brief Get parametric coordinates of element nodes @@ -413,11 +438,11 @@ struct LinearTetrahedronShape { */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector3 const &xi) const { - assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); + assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1] - xi[2]; - assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); - assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); // clang-format off return {L0, xi[0], //L1 @@ -429,11 +454,12 @@ struct LinearTetrahedronShape { /** * @brief Get gradients of shape functions in parametric coordinates * @detail gradients are constant - ignoring parametric coord input - * @return Array of constant gradients [dN0/dxi0, dN0/dxi1, dN0/dxi2, dN1/dxi0, ...] - * Constant gradients: [[-1,-1,-1], [1,0,0], [0,1,0], [0,0,1]] + * @return Array of constant gradients [dN0/dxi0, dN0/dxi1, dN0/dxi2, + * dN1/dxi0, ...] Constant gradients: [[-1,-1,-1], [1,0,0], [0,1,0], [0,0,1]] */ KOKKOS_INLINE_FUNCTION - Kokkos::Array getLocalGradients(Vector3 const &) const { + Kokkos::Array + getLocalGradients(Vector3 const &) const { // clang-format off return {-1, -1, -1, 1, 0, 0, @@ -459,16 +485,19 @@ struct LinearTetrahedronShape { * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct QuadraticTetrahedronShape { - static const size_t numNodes = 10; ///< Number of nodes (10 for quadratic tetrahedron) - 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] = {4, 6}; ///< 4 vertices, 6 edges - constexpr static size_t DofsPerHolder[2] = {1, 1}; ///< 1 DOF per vertex/edge - constexpr static size_t Order = 2; ///< Polynomial order (quadratic) + static const size_t numNodes = + 10; ///< Number of nodes (10 for quadratic tetrahedron) + 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] = {4, 6}; ///< 4 vertices, 6 edges + constexpr static size_t DofsPerHolder[2] = {1, 1}; ///< 1 DOF per vertex/edge + constexpr static size_t Order = 2; ///< Polynomial order (quadratic) /** * @brief Get parametric coordinates of element nodes - * @return Array of node coordinates for 10 nodes (4 at vertices, 6 at edge midpoints) + * @return Array of node coordinates for 10 nodes (4 at vertices, 6 at edge + * midpoints) */ KOKKOS_INLINE_FUNCTION Kokkos::Array getNodeParametricCoords() const { @@ -496,15 +525,16 @@ struct QuadraticTetrahedronShape { * @return Array of 10 shape function values for quadratic tetrahedron * N_i = L_i(2L_i - 1) for vertices (i=0,1,2,3) * N_i = 4*L_j*L_k for edge midpoints (i=4-9) - * where L0=1-xi0-xi1-xi2, L1=xi0, L2=xi1, L3=xi2 are barycentric coordinates + * where L0=1-xi0-xi1-xi2, L1=xi0, L2=xi1, L3=xi2 are barycentric + * coordinates */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector3 const &xi) const { - assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); + assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1] - xi[2]; - assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); - assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); const Real L1 = xi[0]; const Real L2 = xi[1]; const Real L3 = xi[2]; @@ -525,16 +555,17 @@ struct QuadraticTetrahedronShape { /** * @brief Get gradients of quadratic shape functions in parametric coordinates * @param xi Parametric coordinates (xi0, xi1, xi2) where xi_i \f$\in\f$ [0,1] - * @return Array of gradient vectors [dN0/dxi0, dN0/dxi1, dN0/dxi2, ..., dN9/dxi0, dN9/dxi1, dN9/dxi2] + * @return Array of gradient vectors [dN0/dxi0, dN0/dxi1, dN0/dxi2, ..., + * dN9/dxi0, dN9/dxi1, dN9/dxi2] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getLocalGradients(Vector3 const &xi) const { - assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); + assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1] - xi[2]; - assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); - assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); + assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); + assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); const Real L1 = xi[0]; const Real L2 = xi[1]; const Real L3 = xi[2]; @@ -560,11 +591,10 @@ struct QuadraticTetrahedronShape { /** * @brief Linear shape functions for 3D wedge (triangular prism) elements * - * 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] + * 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 (node a+3 is joined to node a by an edge along \zeta): * - Node 0 at (0, 0 , -1) @@ -574,15 +604,16 @@ struct QuadraticTetrahedronShape { * - Node 4 at (1, 0 , 1) * - Node 5 at (0, 1 , 1) * - * Shape functions use barycentric coordinates L0 = 1-xi0-xi1, L1 = xi0, L2 = xi1 of the triangle and, for 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$, with - * \f$\zeta_a = -1\f$ for nodes 0-2 and \f$+1\f$ for nodes 3-5 + * 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 @@ -612,7 +643,7 @@ 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] where + * @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 */ diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 90ab0ea..7841e27 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] @@ -276,13 +276,14 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { // clang-format on QuadraticLagrangeShapeField qlsf(meshInfo, {vtxField, edgeField}); return FieldWithController{kk_ctrl, - qlsf}; + qlsf}; //! [CreateFieldControllerQuadratic] } else { fail("CreateLagrangeField does not support the specified " "combination of order %d and dimension %d.\n", order, dim); - return FieldWithController{}; // silence compiler warning + return FieldWithController{}; // silence compiler warning } }; @@ -309,7 +310,8 @@ template auto CreateCoordinateField(const MeshInfo &meshInfo) { // current coordinate fields are assumed to be linear over the elements - return CreateLagrangeField(meshInfo); + return CreateLagrangeField( + meshInfo); } } // namespace MeshField From c738b014ff9fdd22fcf3005d6b12487a9272e05a Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 01:38:16 -0400 Subject: [PATCH 33/45] fix zeta order initialization --- src/MeshField_Integrate.hpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index 43fb0fd..87ce7ef 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -286,14 +286,14 @@ template auto const getIntegration() { } //! [getIntegration] -// wedges use triOrder in (xi0, xi1) and zetaOrder in zeta; other topologies -// use triOrder +// zetaOrder is used only for wedges and defaults to order template -auto getIntegrationPoints(int triOrder, int zetaOrder) { +auto getIntegrationPoints(int order, int zetaOrder = -1) { if constexpr (topo == Wedge) { - return getWedgeIntegrationPoints(triOrder, zetaOrder); + return getWedgeIntegrationPoints(order, zetaOrder < 0 ? order : zetaOrder); } else { - return getIntegrationPoints(triOrder); + auto ip = getIntegration()->getAccurate(order)->getPoints(); + return ip; } } @@ -375,11 +375,11 @@ auto getJacobianDeterminants(FieldElement &fes, class Integrator { public: /** \brief Construct an Integrator given an order of accuracy. */ - Integrator(int o) : order(o) {} + Integrator(int o) : order(o), zetaOrder(o) {} /** \brief Construct an Integrator for wedges, with separate orders of * accuracy for the triangle (xi0, xi1) and zeta directions. */ - Integrator(int triOrder, int zetaOrderIn) - : order(triOrder), zetaOrder(zetaOrderIn) {} + Integrator(int triOrderIn, int zetaOrderIn) + : order(triOrderIn), zetaOrder(zetaOrderIn) {} virtual ~Integrator(){}; /** \brief User callback: process entry 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 34/45] 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 7ef97e15b8f06a820232007fc53c2b80efde2139 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 03:01:16 -0400 Subject: [PATCH 35/45] Clean up the triangle N4 reference comment --- src/MeshField_Integrate.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index c2abdec..ceee260 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -82,6 +82,7 @@ class TriangleIntegration : public EntityIntegration<2> { // copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). // weights sum to one, so divided by 2, the area of the reference // triangle. The p = 3 rule has a negative weight, so order 3 uses N4. + class N4 : public Integration<2> { public: virtual int countPoints() const { return 6; } @@ -101,7 +102,7 @@ class TriangleIntegration : public EntityIntegration<2> { } virtual int getAccuracy() const { return 4; } }; // end N4 - + virtual int countIntegrations() const { return 3; } virtual Integration<2> const *getIntegration(int i) const { static N1 i1; From b0c30177024ef4e4fe06583855ac0ca86b0556df Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 03:06:13 -0400 Subject: [PATCH 36/45] Revert "clang-format the src files" This reverts commit 33376df7dcf24a90f7106f41c68e1c421f6a9436. --- src/CabanaController.hpp | 4 +- src/MeshField.hpp | 4 +- src/MeshField_Config.hpp.in | 2 +- src/MeshField_Element.hpp | 10 +- src/MeshField_Shape.hpp | 213 +++++++++++++++-------------------- src/MeshField_ShapeField.hpp | 14 +-- 6 files changed, 106 insertions(+), 141 deletions(-) diff --git a/src/CabanaController.hpp b/src/CabanaController.hpp index 289ce53..aadf534 100644 --- a/src/CabanaController.hpp +++ b/src/CabanaController.hpp @@ -111,8 +111,8 @@ class CabanaController { template using wrapper_slice_t = CabanaSliceWrapper, T>; - using AoSoATuple = std::tuple< - Cabana::AoSoA, MemorySpace, vecLen>...>; + using AoSoATuple = + std::tuple, MemorySpace, vecLen>...>; template static auto construct(std::integer_sequence, diff --git a/src/MeshField.hpp b/src/MeshField.hpp index a4bd8f6..dd541b8 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -223,8 +223,8 @@ struct QuadraticTetrahedronToField { const auto tetDim = 3; const auto vtxDim = 0; const auto ignored = -1; - // cyclic rotation of the omegah vertex order to map to the meshfields - // order defined by the shape functions in MeshField_Shape.hpp + // cyclic rotation of the omegah vertex order to map to the meshfields order + // defined by the shape functions in MeshField_Shape.hpp const auto localVtxIdx = (Omega_h::simplex_down_template( tetDim, vtxDim, dofHolderIdx, ignored) + 3) % diff --git a/src/MeshField_Config.hpp.in b/src/MeshField_Config.hpp.in index 5287cd1..cb05f97 100644 --- a/src/MeshField_Config.hpp.in +++ b/src/MeshField_Config.hpp.in @@ -1,2 +1,2 @@ #cmakedefine MeshFields_USE_EXCEPTIONS -#cmakedefine MeshFields_VERSION @MeshFields_VERSION @ +#cmakedefine MeshFields_VERSION @MeshFields_VERSION@ diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index 6ed5b1b..cb2890d 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -88,19 +88,17 @@ KOKKOS_INLINE_FUNCTION auto addTensorProduct(VecA const &a, VecB const &b, namespace MeshField { + // /** * @brief Supports mapping between mesh (i.e., Omega_h) ordering and MeshFields * ordering * - * Return type used by structs/classes that implement the - * ElementDofHolderAccessor parenthesis operator. See @ref - * adding-shape-functions-omegah "Adding Shape Functions Omega_h Element - * Topologies" + * Return type used by structs/classes that implement the ElementDofHolderAccessor parenthesis operator. See @ref adding-shape-functions-omegah "Adding Shape Functions Omega_h Element Topologies" */ struct ElementToDofHolderMap { - /** Which node is being accessed; non-zero + /** Which node is being accessed; non-zero * for mesh entities associated with multiple nodes */ LO node; @@ -116,7 +114,7 @@ struct ElementToDofHolderMap { */ LO entity; - /** Topological type of the entity at entityIndex + /** Topological type of the entity at entityIndex * (e.g., Vertex, Edge, Triangle) */ Mesh_Topology topo; diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index 0b0195b..8489814 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -23,11 +23,8 @@ namespace { * @param val Reference value * @return true if xi >= val within machine precision */ -KOKKOS_INLINE_FUNCTION bool greaterThanOrEqual(MeshField::Real xi, - const MeshField::Real val, - const MeshField::Real tol) { - if (xi > val) - return true; +KOKKOS_INLINE_FUNCTION bool greaterThanOrEqual(MeshField::Real xi, const MeshField::Real val, const MeshField::Real tol) { + if ( xi > val ) return true; return (Kokkos::fabs(xi - val) <= tol); } @@ -39,12 +36,10 @@ KOKKOS_INLINE_FUNCTION bool greaterThanOrEqual(MeshField::Real xi, * @return true if all xi[i] >= val within machine precision */ template -KOKKOS_INLINE_FUNCTION bool eachGreaterThanOrEqual(Array &xi, - const MeshField::Real val, - const MeshField::Real tol) { +KOKKOS_INLINE_FUNCTION bool eachGreaterThanOrEqual(Array &xi, const MeshField::Real val, const MeshField::Real tol) { auto gt = true; for (size_t i = 0; i < xi.size(); i++) { - gt = gt && greaterThanOrEqual(xi[i], val, tol); + gt = gt && greaterThanOrEqual(xi[i],val,tol); } return gt; } @@ -55,11 +50,8 @@ KOKKOS_INLINE_FUNCTION bool eachGreaterThanOrEqual(Array &xi, * @param val Reference value * @return true if xi <= val within machine precision */ -KOKKOS_INLINE_FUNCTION bool lessThanOrEqual(MeshField::Real xi, - const MeshField::Real val, - const MeshField::Real tol) { - if (xi < val) - return true; +KOKKOS_INLINE_FUNCTION bool lessThanOrEqual(MeshField::Real xi, const MeshField::Real val, const MeshField::Real tol) { + if ( xi < val ) return true; return (Kokkos::fabs(xi - val) <= tol); } @@ -71,12 +63,10 @@ KOKKOS_INLINE_FUNCTION bool lessThanOrEqual(MeshField::Real xi, * @return true if all xi[i] <= val within machine precision */ template -KOKKOS_INLINE_FUNCTION bool eachLessThanOrEqual(Array &xi, - const MeshField::Real val, - const MeshField::Real tol) { +KOKKOS_INLINE_FUNCTION bool eachLessThanOrEqual(Array &xi, const MeshField::Real val, const MeshField::Real tol) { auto lt = true; for (size_t i = 0; i < xi.size(); i++) { - lt = lt && lessThanOrEqual(xi[i], val, tol); + lt = lt && lessThanOrEqual(xi[i],val,tol); } return lt; } @@ -107,11 +97,10 @@ using Vector3 = Kokkos::Array; * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct LinearEdgeShape { - static const size_t numNodes = 2; ///< Number of nodes (2 for linear edge) - static const size_t meshEntDim = 1; ///< Mesh entity dimension (1D) - constexpr static Mesh_Topology DofHolders[1] = { - Vertex}; ///< DOFs located at vertices - constexpr static size_t Order = 1; ///< Polynomial order (linear) + static const size_t numNodes = 2; ///< Number of nodes (2 for linear edge) + static const size_t meshEntDim = 1; ///< Mesh entity dimension (1D) + constexpr static Mesh_Topology DofHolders[1] = {Vertex}; ///< DOFs located at vertices + constexpr static size_t Order = 1; ///< Polynomial order (linear) /** * @brief Get parametric coordinates of element nodes @@ -132,8 +121,8 @@ struct LinearEdgeShape { */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector1 const &xi) const { - assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi, -1.0, ParametricCoordTol)); + assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi,-1.0,ParametricCoordTol)); // clang-format off return {(1.0 - xi[0]) / 2.0, (1.0 + xi[0]) / 2.0}; @@ -165,26 +154,22 @@ struct LinearEdgeShape { * - Node 1 at (1, 0) * - Node 2 at (0, 1) * - * Shape functions use barycentric coordinates where L0 = 1-xi0-xi1, L1 = xi0, - * L2 = xi1 + * Shape functions use barycentric coordinates where L0 = 1-xi0-xi1, L1 = xi0, L2 = xi1 * * Shape functions and ordering from: * Zienkiewicz, Taylor, and Zhu * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct LinearTriangleShape { - static const size_t numNodes = 3; ///< Number of nodes (3 for linear triangle) - static const size_t numComponentsPerDof = - 1; ///< Components per DOF (scalar field) - static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) - constexpr static Mesh_Topology DofHolders[1] = { - Vertex}; ///< DOFs located at vertices - constexpr static size_t Order = 1; ///< Polynomial order (linear) + static const size_t numNodes = 3; ///< Number of nodes (3 for linear triangle) + static const size_t numComponentsPerDof = 1; ///< Components per DOF (scalar field) + static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) + constexpr static Mesh_Topology DofHolders[1] = {Vertex}; ///< DOFs located at vertices + constexpr static size_t Order = 1; ///< Polynomial order (linear) /** * @brief Get parametric coordinates of element nodes - * @return Array of node coordinates: [node0_xi0, node0_xi1, node1_xi0, - * node1_xi1, ...] + * @return Array of node coordinates: [node0_xi0, node0_xi1, node1_xi0, node1_xi1, ...] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getNodeParametricCoords() const { @@ -203,11 +188,11 @@ struct LinearTriangleShape { */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector2 const &xi) const { - assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); + assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1]; - assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); - assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); + assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); + assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); // clang-format off return {L0, xi[0], //L1 @@ -218,12 +203,11 @@ struct LinearTriangleShape { /** * @brief Get gradients of shape functions * @detail gradients are constant - ignoring parametric coord input - * @return Array of gradients [dN0/dxi0, dN0/dxi1, dN1/dxi0, dN1/dxi1, - * dN2/dxi0, dN2/dxi1] Constant gradients: [[-1,-1], [1,0], [0,1]] + * @return Array of gradients [dN0/dxi0, dN0/dxi1, dN1/dxi0, dN1/dxi1, dN2/dxi0, dN2/dxi1] + * Constant gradients: [[-1,-1], [1,0], [0,1]] */ KOKKOS_INLINE_FUNCTION - Kokkos::Array - getLocalGradients(Vector2 const &) const { + Kokkos::Array getLocalGradients(Vector2 const &) const { // clang-format off return { -1,-1, //first vector 1, 0, @@ -244,16 +228,14 @@ struct LinearTriangleShape { * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct LinearTriangleCoordinateShape { - static const size_t numNodes = 3; ///< Number of nodes - static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) - constexpr static Mesh_Topology DofHolders[1] = { - Vertex}; ///< DOFs located at vertices - constexpr static size_t Order = 1; ///< Polynomial order (linear) + static const size_t numNodes = 3; ///< Number of nodes + static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) + constexpr static Mesh_Topology DofHolders[1] = {Vertex}; ///< DOFs located at vertices + constexpr static size_t Order = 1; ///< Polynomial order (linear) /** * @brief Get parametric coordinates of element nodes - * @return Array of node coordinates: [node0_xi0, node0_xi1, node1_xi0, - * node1_xi1, node2_xi0, node2_xi1] + * @return Array of node coordinates: [node0_xi0, node0_xi1, node1_xi0, node1_xi1, node2_xi0, node2_xi1] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getNodeParametricCoords() const { @@ -267,16 +249,15 @@ struct LinearTriangleCoordinateShape { /** * @brief Evaluate shape functions at parametric coordinate * @param xi Parametric coordinates (xi0, xi1) where xi0, xi1 \f$\in\f$ [0,1] - * @return Array of shape function values using barycentric coordinates [L0, - * L1, L2] + * @return Array of shape function values using barycentric coordinates [L0, L1, L2] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector2 const &xi) const { - assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); + assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1]; - assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); - assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); + assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); + assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); // clang-format off return {L0, xi[0], @@ -302,19 +283,16 @@ struct LinearTriangleCoordinateShape { */ //! [QuadraticTriangleShape] struct QuadraticTriangleShape { - static const size_t numNodes = - 6; ///< Number of nodes (6 for quadratic triangle) - static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) - constexpr static Mesh_Topology DofHolders[2] = { - Vertex, Edge}; ///< DOFs at vertices and edges - constexpr static size_t NumDofHolders[2] = {3, 3}; ///< 3 vertices, 3 edges - constexpr static size_t DofsPerHolder[2] = {1, 1}; ///< 1 DOF per vertex/edge - constexpr static size_t Order = 2; ///< Polynomial order (quadratic) + static const size_t numNodes = 6; ///< Number of nodes (6 for quadratic triangle) + static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) + constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge}; ///< DOFs at vertices and edges + constexpr static size_t NumDofHolders[2] = {3, 3}; ///< 3 vertices, 3 edges + constexpr static size_t DofsPerHolder[2] = {1, 1}; ///< 1 DOF per vertex/edge + constexpr static size_t Order = 2; ///< Polynomial order (quadratic) /** * @brief Get parametric coordinates of element nodes - * @return Array of node coordinates for 6 nodes (3 at vertices, 3 at edge - * midpoints) + * @return Array of node coordinates for 6 nodes (3 at vertices, 3 at edge midpoints) */ KOKKOS_INLINE_FUNCTION Kokkos::Array getNodeParametricCoords() const { @@ -342,11 +320,11 @@ struct QuadraticTriangleShape { */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector2 const &xi) const { - assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); + assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1]; - assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); - assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); + assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); + assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); const Real L1 = xi[0]; const Real L2 = xi[1]; // clang-format off @@ -362,17 +340,16 @@ struct QuadraticTriangleShape { /** * @brief Get gradients of quadratic shape functions at parametric coordinate * @param xi Parametric coordinates (xi0, xi1) where xi0, xi1 \f$\in\f$ [0,1] - * @return Array of gradient vectors [dN0/dxi0, dN0/dxi1, ..., dN5/dxi0, - * dN5/dxi1] + * @return Array of gradient vectors [dN0/dxi0, dN0/dxi1, ..., dN5/dxi0, dN5/dxi1] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getLocalGradients(Vector2 const &xi) const { - assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); + assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1]; - assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); - assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); + assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); + assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); const Real L1 = xi[0]; const Real L2 = xi[1]; // clang-format off @@ -409,12 +386,10 @@ struct QuadraticTriangleShape { * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct LinearTetrahedronShape { - static const size_t numNodes = - 4; ///< Number of nodes (4 for linear tetrahedron) - 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) + static const size_t numNodes = 4; ///< Number of nodes (4 for linear tetrahedron) + 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) /** * @brief Get parametric coordinates of element nodes @@ -438,11 +413,11 @@ struct LinearTetrahedronShape { */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector3 const &xi) const { - assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); + assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1] - xi[2]; - assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); - assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); + assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); + assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); // clang-format off return {L0, xi[0], //L1 @@ -454,12 +429,11 @@ struct LinearTetrahedronShape { /** * @brief Get gradients of shape functions in parametric coordinates * @detail gradients are constant - ignoring parametric coord input - * @return Array of constant gradients [dN0/dxi0, dN0/dxi1, dN0/dxi2, - * dN1/dxi0, ...] Constant gradients: [[-1,-1,-1], [1,0,0], [0,1,0], [0,0,1]] + * @return Array of constant gradients [dN0/dxi0, dN0/dxi1, dN0/dxi2, dN1/dxi0, ...] + * Constant gradients: [[-1,-1,-1], [1,0,0], [0,1,0], [0,0,1]] */ KOKKOS_INLINE_FUNCTION - Kokkos::Array - getLocalGradients(Vector3 const &) const { + Kokkos::Array getLocalGradients(Vector3 const &) const { // clang-format off return {-1, -1, -1, 1, 0, 0, @@ -485,19 +459,16 @@ struct LinearTetrahedronShape { * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ struct QuadraticTetrahedronShape { - static const size_t numNodes = - 10; ///< Number of nodes (10 for quadratic tetrahedron) - 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] = {4, 6}; ///< 4 vertices, 6 edges - constexpr static size_t DofsPerHolder[2] = {1, 1}; ///< 1 DOF per vertex/edge - constexpr static size_t Order = 2; ///< Polynomial order (quadratic) + static const size_t numNodes = 10; ///< Number of nodes (10 for quadratic tetrahedron) + 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] = {4, 6}; ///< 4 vertices, 6 edges + constexpr static size_t DofsPerHolder[2] = {1, 1}; ///< 1 DOF per vertex/edge + constexpr static size_t Order = 2; ///< Polynomial order (quadratic) /** * @brief Get parametric coordinates of element nodes - * @return Array of node coordinates for 10 nodes (4 at vertices, 6 at edge - * midpoints) + * @return Array of node coordinates for 10 nodes (4 at vertices, 6 at edge midpoints) */ KOKKOS_INLINE_FUNCTION Kokkos::Array getNodeParametricCoords() const { @@ -525,16 +496,15 @@ struct QuadraticTetrahedronShape { * @return Array of 10 shape function values for quadratic tetrahedron * N_i = L_i(2L_i - 1) for vertices (i=0,1,2,3) * N_i = 4*L_j*L_k for edge midpoints (i=4-9) - * where L0=1-xi0-xi1-xi2, L1=xi0, L2=xi1, L3=xi2 are barycentric - * coordinates + * where L0=1-xi0-xi1-xi2, L1=xi0, L2=xi1, L3=xi2 are barycentric coordinates */ KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector3 const &xi) const { - assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); + assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1] - xi[2]; - assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); - assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); + assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); + assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); const Real L1 = xi[0]; const Real L2 = xi[1]; const Real L3 = xi[2]; @@ -555,17 +525,16 @@ struct QuadraticTetrahedronShape { /** * @brief Get gradients of quadratic shape functions in parametric coordinates * @param xi Parametric coordinates (xi0, xi1, xi2) where xi_i \f$\in\f$ [0,1] - * @return Array of gradient vectors [dN0/dxi0, dN0/dxi1, dN0/dxi2, ..., - * dN9/dxi0, dN9/dxi1, dN9/dxi2] + * @return Array of gradient vectors [dN0/dxi0, dN0/dxi1, dN0/dxi2, ..., dN9/dxi0, dN9/dxi1, dN9/dxi2] */ KOKKOS_INLINE_FUNCTION Kokkos::Array getLocalGradients(Vector3 const &xi) const { - assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); - assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); + assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol)); + assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol)); const Real L0 = 1 - xi[0] - xi[1] - xi[2]; - assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); - assert(lessThanOrEqual(L0, 1.0, ParametricCoordTol)); + assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol)); + assert(lessThanOrEqual(L0,1.0,ParametricCoordTol)); const Real L1 = xi[0]; const Real L2 = xi[1]; const Real L3 = xi[2]; @@ -591,10 +560,11 @@ struct QuadraticTetrahedronShape { /** * @brief Linear shape functions for 3D wedge (triangular prism) elements * - * 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] + * 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 (node a+3 is joined to node a by an edge along \zeta): * - Node 0 at (0, 0 , -1) @@ -604,16 +574,15 @@ struct QuadraticTetrahedronShape { * - Node 4 at (1, 0 , 1) * - Node 5 at (0, 1 , 1) * - * Shape functions use barycentric coordinates L0 = 1-xi0-xi1, L1 = xi0, L2 = - * xi1 of the triangle and, for 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$, with - * \f$\zeta_a = -1\f$ for nodes 0-2 and \f$+1\f$ for nodes 3-5 + * 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 @@ -643,7 +612,7 @@ 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] where + * @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 */ diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 7841e27..90ab0ea 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] @@ -276,14 +276,13 @@ auto CreateLagrangeField(const MeshInfo &meshInfo) { // clang-format on QuadraticLagrangeShapeField qlsf(meshInfo, {vtxField, edgeField}); return FieldWithController{kk_ctrl, - qlsf}; + qlsf}; //! [CreateFieldControllerQuadratic] } else { fail("CreateLagrangeField does not support the specified " "combination of order %d and dimension %d.\n", order, dim); - return FieldWithController{}; // silence compiler warning + return FieldWithController{}; // silence compiler warning } }; @@ -310,8 +309,7 @@ template auto CreateCoordinateField(const MeshInfo &meshInfo) { // current coordinate fields are assumed to be linear over the elements - return CreateLagrangeField( - meshInfo); + return CreateLagrangeField(meshInfo); } } // namespace MeshField From 9011685d09a8bb075c9bd96da9fc495a6129d065 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 03:21:59 -0400 Subject: [PATCH 37/45] Keep the triangle and tetrahedron blocks as on triangle-quadrature --- src/MeshField_Integrate.hpp | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index 87ce7ef..1c92f24 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -47,10 +47,11 @@ template class EntityIntegration { virtual int countIntegrations() const = 0; virtual Integration const *getIntegration(int i) const = 0; }; -// Gauss-Legendre rules on [-1, 1] from -// N1 to N6 are copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc -// (EdgeIntegration, lines 27-129). The n-point rule is exact for degree -// 2n - 1, and the weights sum to 2, the length of the reference edge. + +// Gauss-Legendre rules on [-1, 1]. N1 to N6 are copied from SCOREC/core +// @ 7cd76473 apf/apfIntegrate.cc (EdgeIntegration, lines 27-129). The n-point +// rule is exact for degree 2n - 1, and the weights sum to 2, the length of +// the reference edge. //! [EdgeIntegration] class EdgeIntegration : public EntityIntegration<1> { public: @@ -160,12 +161,12 @@ class TriangleIntegration : public EntityIntegration<2> { public: virtual int countPoints() const { return 3; } virtual std::vector> getPoints() const { - return {IntegrationPoint(Vector2{0.166666666666667, 0.166666666666667}, - 1. / 3. / 2.0, 2, 0), - IntegrationPoint(Vector2{0.666666666666667, 0.166666666666667}, - 1. / 3. / 2.0, 2, 0), - IntegrationPoint(Vector2{0.166666666666667, 0.666666666666667}, - 1. / 3. / 2.0, 2, 0)}; + return {IntegrationPoint( + Vector2{0.166666666666667, 0.166666666666667}, 1. / 3. / 2.0, 2, 0), + IntegrationPoint( + Vector2{0.666666666666667, 0.166666666666667}, 1. / 3. / 2.0, 2, 0), + IntegrationPoint( + Vector2{0.166666666666667, 0.666666666666667}, 1. / 3. / 2.0, 2, 0)}; } virtual int getAccuracy() const { return 2; } }; // end N2 @@ -173,10 +174,10 @@ class TriangleIntegration : public EntityIntegration<2> { // Symmetric Gaussian quadrature rules for the triangle from // D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature // rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). - // N4 is p = 4 rules of from the Appendix II (p. 1140), - // copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). - // weights sum to one, so divided by 2, the area of the reference - // triangle. The p = 3 rule has a negative weight, so order 3 uses N4. + // N4 is the p = 4 rule of Appendix II (p. 1140), copied from + // SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). The weights + // sum to one, so they are divided by 2, the area of the reference triangle. + // The p = 3 rule has a negative weight, so order 3 uses N4. class N4 : public Integration<2> { public: virtual int countPoints() const { return 6; } @@ -208,7 +209,6 @@ class TriangleIntegration : public EntityIntegration<2> { }; //! [TriangleIntegration] -//! [TetrahedronIntegrationIntegration] class TetrahedronIntegration : public EntityIntegration<3> { public: class N1 : public Integration<3> { @@ -247,7 +247,6 @@ class TetrahedronIntegration : public EntityIntegration<3> { return integrations[i]; } }; -//! [TetrahedronIntegrationIntegration] // Wedge rules are the tensor product of a triangle rule in (xi0, xi1) and an // edge rule in zeta. The rule for (triOrder, zetaOrder) integrates From 7fc9e565e62c2c34d4f3963934156333fb49bca5 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Tue, 27 Oct 2026 03:33:02 -0400 Subject: [PATCH 38/45] Fix the N4 comment wording and indentation --- src/MeshField_Integrate.hpp | 15 +++++++-------- 1 file changed, 7 insertions(+), 8 deletions(-) diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index ceee260..e370fe0 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -75,14 +75,13 @@ class TriangleIntegration : public EntityIntegration<2> { virtual int getAccuracy() const { return 2; } }; // end N2 -// Symmetric Gaussian quadrature rules for the triangle from -// D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature -// rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). -// N4 is p = 4 rules of from the Appendix II (p. 1140), -// copied from SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). -// weights sum to one, so divided by 2, the area of the reference -// triangle. The p = 3 rule has a negative weight, so order 3 uses N4. - + // Symmetric Gaussian quadrature rules for the triangle from + // D. A. Dunavant, "High degree efficient symmetrical Gaussian quadrature + // rules for the triangle", Int. J. Numer. Meth. Eng. 21, 1129-1148 (1985). + // N4 is the p = 4 rule of Appendix II (p. 1140), copied from + // SCOREC/core @ 7cd76473 apf/apfIntegrate.cc (lines 175-190). The weights + // sum to one, so they are divided by 2, the area of the reference triangle. + // The p = 3 rule has a negative weight, so order 3 uses N4. class N4 : public Integration<2> { public: virtual int countPoints() const { return 6; } 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 39/45] 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 82b6e1b569e47c1d434a4d5c6be445f9949ec59f Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 28 Oct 2026 11:07:42 -0400 Subject: [PATCH 40/45] Adapt TriangleIntegrationTest to per-element integration points --- test/testTriangleIntegration.cpp | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/test/testTriangleIntegration.cpp b/test/testTriangleIntegration.cpp index 5749068..8aacc8d 100644 --- a/test/testTriangleIntegration.cpp +++ b/test/testTriangleIntegration.cpp @@ -76,15 +76,17 @@ class MonomialIntegrator : public MeshField::Integrator { : MeshField::Integrator(order), fe(feIn), px(pxIn), py(pyIn) {} void atPoints(Kokkos::View p, Kokkos::View w, Kokkos::View dV) override { - const int numPtsPerElem = p.extent(0) / fe.numMeshEnts; - const auto x = MeshField::evaluate(fe, p, numPtsPerElem); + // p and w hold one entry per point of a single element; x and dV hold + // one entry per point of every element + const auto x = MeshField::evaluateStaticLength(fe, p); + const int numPtsPerElem = w.extent(0); const Real a = px, b = py; Real sum = 0; Kokkos::parallel_reduce( - "integrateMonomial", w.extent(0), + "integrateMonomial", dV.extent(0), KOKKOS_LAMBDA(const int i, Real &lsum) { - lsum += - w(i) * dV(i) * Kokkos::pow(x(i, 0), a) * Kokkos::pow(x(i, 1), b); + lsum += w(i % numPtsPerElem) * dV(i) * Kokkos::pow(x(i, 0), a) * + Kokkos::pow(x(i, 1), b); }, sum); integral += sum; 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 41/45] 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)) { From 43556ed0ae1ab5b30c658d3d4ac0efc82c701d4d Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 28 Oct 2026 11:11:45 -0400 Subject: [PATCH 42/45] Adapt WedgeIntegrationTest to per-element integration points --- test/testWedgeIntegration.cpp | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/test/testWedgeIntegration.cpp b/test/testWedgeIntegration.cpp index 65b36d2..063175d 100644 --- a/test/testWedgeIntegration.cpp +++ b/test/testWedgeIntegration.cpp @@ -103,10 +103,16 @@ class VolumeIntegrator : public MeshField::Integrator { : MeshField::Integrator(triOrder, zetaOrder) {} void atPoints(Kokkos::View, Kokkos::View w, Kokkos::View dV) override { + // w holds one weight per point of a single element; dV holds one entry + // per point of every element + const int numPtsPerElem = w.extent(0); Real sum = 0; Kokkos::parallel_reduce( - "volume", w.extent(0), - KOKKOS_LAMBDA(const int i, Real &s) { s += w(i) * dV(i); }, sum); + "volume", dV.extent(0), + KOKKOS_LAMBDA(const int i, Real &s) { + s += w(i % numPtsPerElem) * dV(i); + }, + sum); volume = sum; } Real volume = 0; @@ -136,8 +142,8 @@ template 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"); } int checkIntegrator() { From adefeb34074f2c36918ca84c53deed951402232f Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 28 Oct 2026 16:05:35 -0400 Subject: [PATCH 43/45] Add extruded wedge element map --- src/MeshField.hpp | 183 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 183 insertions(+) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index a967837..9f05e1b 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -284,6 +284,189 @@ template auto getTetrahedronElement(Omega_h::Mesh &mesh) { } } +/** + * @brief Maps the nodes of LinearWedgeShape wedges, extruded from a 2D + * triangle mesh, to vertices + * @details + * The 2D mesh is stacked into numLayers layers and every plane has the same + * triangulation. Vertex v on plane p is vertex p * nVtx2d + v. Wedge w is in + * layer k = w / nTri2d and comes from triangle t = w % nTri2d. Layer k runs + * from plane k to plane k + 1, or to plane (k + 1) % numLayers when periodic. + * Wedge corner a is local vertex a of triangle t, or local vertex {0, 2, 1}[a] + * when reverseOrientation is true. Node a is corner a on plane k and node + * a + 3 is corner a on the top plane of the layer. + * + * With coordinates (R, phi, Z), phi increasing with the layer and triangles + * counterclockwise in (R, Z), det J < 0. Setting reverseOrientation = true + * makes it positive. + */ +struct ExtrudedLinearWedgeToVertexField { + Omega_h::LOs triVerts; + int nVtx2d; + int nTri2d; + int numLayers; + bool periodic; + bool reverseOrientation; + ExtrudedLinearWedgeToVertexField(Omega_h::Mesh &mesh2d, int numLayersIn, + bool periodicIn, bool reverseOrientationIn) + : triVerts(mesh2d.ask_elem_verts()), nVtx2d(mesh2d.nverts()), + nTri2d(mesh2d.nfaces()), numLayers(numLayersIn), periodic(periodicIn), + reverseOrientation(reverseOrientationIn) { + if (mesh2d.dim() != 2 || mesh2d.family() != OMEGA_H_SIMPLEX) { + MeshField::fail( + "The mesh passed to %s must be 2D and simplex (triangles)\n", + __func__); + } + if (numLayers < 1 || (periodic && numLayers < 2)) { + MeshField::fail("%s needs numLayers >= 1, and >= 2 when periodic; got " + "%d\n", + __func__, numLayers); + } + } + + 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 layer = wedge / nTri2d; + const MeshField::LO tri = wedge % nTri2d; + const MeshField::LO corner = wedgeNodeIdx % 3; + // reverseOrientation swaps corners 1 and 2 + const MeshField::LO triNodeIdx = + (reverseOrientation && corner != 0) ? 3 - corner : corner; + const auto triDim = 2; + const auto vtxDim = 0; + const auto ignored = -1; + const auto localVtxIdx = + Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored); + const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim); + const MeshField::LO vtx2d = triVerts[(tri * triToVtxDegree) + localVtxIdx]; + const MeshField::LO top = periodic ? (layer + 1) % numLayers : layer + 1; + const MeshField::LO plane = (wedgeNodeIdx < 3) ? layer : top; + return {0, wedgeCompIdx, plane * nVtx2d + vtx2d, MeshField::Vertex}; + } +}; + +/** + * @brief Maps the nodes of LinearCubicWedgeShape wedges, extruded from a 2D + * triangle mesh, to vertices and vertical edges + * @details + * Nodes 0-5 map as in ExtrudedLinearWedgeToVertexField. The vertical edge + * above 2D vertex v in layer k is edge k * nVtx2d + v. Nodes 6 + 2a and + * 7 + 2a are DOFs 0 (zeta = -1/3) and 1 (zeta = +1/3) of the edge above + * corner a. The orientation note of ExtrudedLinearWedgeToVertexField applies. + */ +struct ExtrudedLinearCubicWedgeToField { + Omega_h::LOs triVerts; + int nVtx2d; + int nTri2d; + int numLayers; + bool periodic; + bool reverseOrientation; + ExtrudedLinearCubicWedgeToField(Omega_h::Mesh &mesh2d, int numLayersIn, + bool periodicIn, bool reverseOrientationIn) + : triVerts(mesh2d.ask_elem_verts()), nVtx2d(mesh2d.nverts()), + nTri2d(mesh2d.nfaces()), numLayers(numLayersIn), periodic(periodicIn), + reverseOrientation(reverseOrientationIn) { + if (mesh2d.dim() != 2 || mesh2d.family() != OMEGA_H_SIMPLEX) { + MeshField::fail( + "The mesh passed to %s must be 2D and simplex (triangles)\n", + __func__); + } + if (numLayers < 1 || (periodic && numLayers < 2)) { + MeshField::fail("%s needs numLayers >= 1, and >= 2 when periodic; got " + "%d\n", + __func__, numLayers); + } + } + + 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 layer = wedge / nTri2d; + const MeshField::LO tri = wedge % nTri2d; + const bool isVtx = wedgeNodeIdx < 6; + const MeshField::LO corner = + isVtx ? wedgeNodeIdx % 3 : (wedgeNodeIdx - 6) / 2; + // reverseOrientation swaps corners 1 and 2 + const MeshField::LO triNodeIdx = + (reverseOrientation && corner != 0) ? 3 - corner : corner; + const auto triDim = 2; + const auto vtxDim = 0; + const auto ignored = -1; + const auto localVtxIdx = + Omega_h::simplex_down_template(triDim, vtxDim, triNodeIdx, ignored); + const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim); + const MeshField::LO vtx2d = triVerts[(tri * triToVtxDegree) + localVtxIdx]; + if (isVtx) { + const MeshField::LO top = periodic ? (layer + 1) % numLayers : layer + 1; + const MeshField::LO plane = (wedgeNodeIdx < 3) ? layer : top; + return {0, wedgeCompIdx, plane * nVtx2d + vtx2d, MeshField::Vertex}; + } + return {(wedgeNodeIdx - 6) % 2, wedgeCompIdx, layer * nVtx2d + vtx2d, + MeshField::Edge}; + } +}; + +/** + * @brief MeshInfo of the wedges from getExtrudedWedgeElement + * @details numEdge counts only the vertical edges, which hold the edge DOFs + * of LinearCubicWedgeShape + */ +inline MeshField::MeshInfo getExtrudedWedgeMeshInfo(Omega_h::Mesh &mesh2d, + int numLayers, + bool periodic = false) { + const int numPlanes = periodic ? numLayers : numLayers + 1; + MeshField::MeshInfo meshInfo; + meshInfo.dim = 3; + meshInfo.numVtx = numPlanes * mesh2d.nverts(); + meshInfo.numEdge = numLayers * mesh2d.nverts(); + meshInfo.numWedge = numLayers * mesh2d.nfaces(); + return meshInfo; +} + +/** + * @brief Shape and element map for wedges extruded from a 2D triangle mesh + * @tparam ShapeOrder 1 for LinearWedgeShape, 3 for LinearCubicWedgeShape + * @details See ExtrudedLinearWedgeToVertexField for the numbering and + * orientation. The caller fills the coordinate field, sized by + * getExtrudedWedgeMeshInfo. + */ +template +auto getExtrudedWedgeElement(Omega_h::Mesh &mesh2d, int numLayers, + bool periodic = false, + bool reverseOrientation = false) { + static_assert(ShapeOrder == 1 || ShapeOrder == 3); + if constexpr (ShapeOrder == 1) { + struct result { + MeshField::LinearWedgeShape shp; + ExtrudedLinearWedgeToVertexField map; + }; + return result{MeshField::LinearWedgeShape(), + ExtrudedLinearWedgeToVertexField(mesh2d, numLayers, periodic, + reverseOrientation)}; + } else if constexpr (ShapeOrder == 3) { + struct result { + MeshField::LinearCubicWedgeShape shp; + ExtrudedLinearCubicWedgeToField map; + }; + return result{MeshField::LinearCubicWedgeShape(), + ExtrudedLinearCubicWedgeToField(mesh2d, numLayers, periodic, + reverseOrientation)}; + } +} + } // namespace Omegah /** From 93ffa65a8348eebdafa5470d4e9f93120e0d62d3 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 28 Oct 2026 16:05:35 -0400 Subject: [PATCH 44/45] Add extruded wedge test --- CMakeLists.txt | 2 + test/testExtrudedWedge.cpp | 413 +++++++++++++++++++++++++++++++++++++ 2 files changed, 415 insertions(+) create mode 100644 test/testExtrudedWedge.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 1714ed1..3e16c15 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -226,6 +226,7 @@ meshfields_add_exe(LinearWedgeShapeTest test/testLinearWedgeShape.cpp) meshfields_add_exe(LinearWedgeElementTest test/testLinearWedgeElement.cpp) meshfields_add_exe(TriangleIntegrationTest test/testTriangleIntegration.cpp) meshfields_add_exe(WedgeIntegrationTest test/testWedgeIntegration.cpp) +meshfields_add_exe(ExtrudedWedgeTest test/testExtrudedWedge.cpp) meshfields_add_exe(Interpolate test/testInterpolate.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) @@ -258,6 +259,7 @@ test_func(LinearWedgeShapeTest ./LinearWedgeShapeTest) test_func(LinearWedgeElementTest ./LinearWedgeElementTest) test_func(TriangleIntegrationTest ./TriangleIntegrationTest) test_func(WedgeIntegrationTest ./WedgeIntegrationTest) +test_func(ExtrudedWedgeTest ./ExtrudedWedgeTest) test_func(Interpolate ./Interpolate) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error diff --git a/test/testExtrudedWedge.cpp b/test/testExtrudedWedge.cpp new file mode 100644 index 0000000..d577415 --- /dev/null +++ b/test/testExtrudedWedge.cpp @@ -0,0 +1,413 @@ +#include "KokkosController.hpp" +#include "MeshField.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_For.hpp" +#include "MeshField_Integrate.hpp" +#include "MeshField_ShapeField.hpp" +#include "Omega_h_build.hpp" +#include "Omega_h_library.hpp" +#include +#include +#include +#include +#include + +// Checks the extruded wedge element maps, getExtrudedWedgeElement and +// getExtrudedWedgeMeshInfo, on a 2D triangle mesh stacked into numLayers +// layers: +// 1. the DOF holder of every wedge node, for open and periodic stacks, linear +// and cubic, with and without reverseOrientation, and the MeshInfo counts, +// 2. the field positions at the nodes and interior points of every wedge, +// 3. the volume from the Integrator, +// 4. the sign of det J with (R, phi, Z) = (x, s, y). +// The coordinates come from analytic maps X(x, y, s) that the elements +// reproduce exactly. Checks 2-4 use open stacks: in a periodic stack the last +// layer closes onto plane 0, which a single map X cannot describe. Failures +// are counted explicitly so the test also fails in builds where assert() is +// disabled. + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::LO; +using MeshField::Real; +using Vec3 = std::array; + +namespace { + +constexpr int numLayers = 3; +constexpr Real H = 1.5; // height of the stack +constexpr Real h = H / numLayers; // height of one layer +constexpr Real A = 1.0; // area of the 2D mesh + +bool isClose(Real a, Real b, Real tol) { + return std::fabs(a - b) <= tol * (1 + std::fabs(b)); +} + +// Analytic maps X(x, y, s) of the point (x, y) of the 2D mesh at height s +using Map = Vec3 (*)(Real, Real, Real); +Vec3 straightMap(Real x, Real y, Real s) { return {x, y, s}; } +Vec3 taperedMap(Real x, Real y, Real s) { return {x * (1 + 0.1 * s), y, s}; } +Vec3 curvedMap(Real x, Real y, Real s) { + return {x * (1 + 0.1 * s) + 0.02 * s * s * s, y, s + 0.05 * s * s * s}; +} +// R = x, phi = s, Z = y +Vec3 cylindricalMap(Real x, Real y, Real s) { return {x, s, y}; } + +// Host copies of the 2D mesh +struct HostMesh2d { + HostMesh2d(Omega_h::Mesh &mesh2d) + : triVerts(mesh2d.ask_elem_verts()), coords(mesh2d.coords()), + nVtx2d(mesh2d.nverts()), nTri2d(mesh2d.nfaces()) {} + Omega_h::HostRead triVerts; + Omega_h::HostRead coords; + int nVtx2d; + int nTri2d; + + // the 2D vertex of wedge corner a above triangle t + LO cornerVtx(LO t, int a, bool reverse) const { + const int corner[2][3] = {{0, 1, 2}, {0, 2, 1}}; + return triVerts[t * 3 + corner[reverse][a]]; + } + + Real signedArea(LO t) const { + Real x[3], y[3]; + for (int a = 0; a < 3; ++a) { + const LO v = triVerts[t * 3 + a]; + x[a] = coords[2 * v]; + y[a] = coords[2 * v + 1]; + } + return 0.5 * + ((x[1] - x[0]) * (y[2] - y[0]) - (x[2] - x[0]) * (y[1] - y[0])); + } + + // X at the local point xi of wedge w in an open stack + Vec3 position(Map X, LO w, Vec3 const &xi, bool reverse) const { + const LO k = w / nTri2d, t = w % nTri2d; + const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; + Real x = 0, y = 0; + for (int a = 0; a < 3; ++a) { + const LO v = cornerVtx(t, a, reverse); + x += L[a] * coords[2 * v]; + y += L[a] * coords[2 * v + 1]; + } + return X(x, y, k * h + h * (1 + xi[2]) / 2); + } +}; + +template +auto createCoordField(const MeshField::MeshInfo &info) { + if constexpr (ShapeOrder == 1) { + return MeshField::CreateCoordinateField( + info); + } else { + return MeshField::CreateLagrangeField< + ExecutionSpace, MeshField::KokkosController, Real, 3, 3, 3, + MeshField::LinearCubicWedgeShape>(info); + } +} + +// Writes X into the coordinate field of an open stack: vertex (p, v) at +// X(x_v, y_v, s_p) and, for the cubic field, DOFs 0 and 1 of edge (k, v) at +// X(x_v, y_v, s_k + h/3) and X(x_v, y_v, s_k + 2h/3). +template +void setExtrudedCoords(CoordField field, const HostMesh2d &m, Map X) { + const int numVtx = (numLayers + 1) * m.nVtx2d; + const int numEdges = numLayers * m.nVtx2d; + Kokkos::View vtx_h("vtx_h", numVtx); + Kokkos::View edge_h("edge_h", numEdges); + for (int v = 0; v < m.nVtx2d; ++v) { + const Real x = m.coords[2 * v], y = m.coords[2 * v + 1]; + for (int p = 0; p <= numLayers; ++p) { + const Vec3 pos = X(x, y, p * h); + for (int d = 0; d < 3; ++d) + vtx_h(p * m.nVtx2d + v, d) = pos[d]; + } + for (int k = 0; k < numLayers; ++k) + for (int dof = 0; dof < 2; ++dof) { + const Vec3 pos = X(x, y, k * h + (dof + 1) * h / 3); + for (int d = 0; d < 3; ++d) + edge_h(k * m.nVtx2d + v, dof, d) = pos[d]; + } + } + auto vtx = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), vtx_h); + auto setVtxCoords = KOKKOS_LAMBDA(const int &i) { + for (int d = 0; d < 3; ++d) + field(i, 0, d, MeshField::Vertex) = vtx(i, d); + }; + MeshField::parallel_for( + {0}, {numVtx}, setVtxCoords, "setVtxCoords"); + if constexpr (ShapeOrder == 3) { + auto edge = Kokkos::create_mirror_view_and_copy(ExecutionSpace(), edge_h); + auto setEdgeCoords = 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); + }; + MeshField::parallel_for( + {0}, {numEdges}, setEdgeCoords, "setEdgeCoords"); + } +} + +const std::vector interiorPts = { + {1. / 3, 1. / 3, 0}, {0.2, 0.3, -0.4}, {0.6, 0.1, 0.7}, {0.05, 0.9, 0.95}}; + +// The nodes of the shape, then interior points +template std::vector nodesAndInteriorPoints() { + const auto xi = Shape().getNodeParametricCoords(); + std::vector pts; + for (size_t n = 0; n < Shape::numNodes; ++n) + pts.push_back({xi[3 * n], xi[3 * n + 1], xi[3 * n + 2]}); + pts.insert(pts.end(), interiorPts.begin(), interiorPts.end()); + return pts; +} + +// the same local points in every wedge, stored element-major +Kokkos::View makeLocalCoords(int numWedges, + const std::vector &pts) { + const int numPts = static_cast(pts.size()); + 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); + return lc; +} + +const char *shapeName(int shapeOrder) { + return shapeOrder == 1 ? "linear" : "cubic"; +} + +// Gets the DOF holder of every (wedge, node) from the map on the device and +// compares it with the layout computed on the host. Also checks that the ids +// are in range, that the top of layer k is the bottom of layer k + 1 (and the +// top of the last layer is the bottom of layer 0 when periodic), and the +// MeshInfo counts. +template +int checkIds(Omega_h::Mesh &mesh2d, const HostMesh2d &m, bool periodic, + bool reverse) { + const auto elem = MeshField::Omegah::getExtrudedWedgeElement( + mesh2d, numLayers, periodic, reverse); + const auto map = elem.map; + const int numNodes = decltype(elem.shp)::numNodes; + const int numWedges = numLayers * m.nTri2d; + // node, entity, topology, and 1 if the component is passed through + Kokkos::View holders("holders", numWedges * numNodes, 4); + Kokkos::parallel_for( + "getHolders", + Kokkos::RangePolicy(0, numWedges * numNodes), + KOKKOS_LAMBDA(const int i) { + const LO w = i / numNodes, n = i % numNodes; + const auto holder = map(n, 0, w, MeshField::Wedge); + LO compOk = holder.component == 0; + for (LO c = 1; c < 3; ++c) + compOk = compOk && map(n, c, w, MeshField::Wedge).component == c; + holders(i, 0) = holder.node; + holders(i, 1) = holder.entity; + holders(i, 2) = holder.topo; + holders(i, 3) = compOk; + }); + const auto holders_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), holders); + + const auto info = + MeshField::Omegah::getExtrudedWedgeMeshInfo(mesh2d, numLayers, periodic); + const int numPlanes = periodic ? numLayers : numLayers + 1; + int idFailed = 0, rangeFailed = 0, stackFailed = 0, countFailed = 0; + for (LO w = 0; w < numWedges; ++w) { + const LO k = w / m.nTri2d, t = w % m.nTri2d; + const LO top = periodic ? (k + 1) % numLayers : k + 1; + for (int n = 0; n < numNodes; ++n) { + const int i = w * numNodes + n; + const LO dof = holders_h(i, 0), ent = holders_h(i, 1), + topo = holders_h(i, 2); + LO expDof = 0, expEnt = 0, expTopo = MeshField::Vertex; + if (n < 6) { + const LO plane = n < 3 ? k : top; + expEnt = plane * m.nVtx2d + m.cornerVtx(t, n % 3, reverse); + } else { + expDof = (n - 6) % 2; + expEnt = k * m.nVtx2d + m.cornerVtx(t, (n - 6) / 2, reverse); + expTopo = MeshField::Edge; + } + idFailed += !(dof == expDof && ent == expEnt && topo == expTopo && + holders_h(i, 3) == 1); + if (topo == MeshField::Vertex) + rangeFailed += !(dof == 0 && ent >= 0 && ent < info.numVtx); + else + rangeFailed += !(topo == MeshField::Edge && (dof == 0 || dof == 1) && + ent >= 0 && ent < info.numEdge); + } + } + for (LO k = 0; k < numLayers; ++k) { + if (k + 1 == numLayers && !periodic) + continue; + const LO next = (k + 1) % numLayers; + for (LO t = 0; t < m.nTri2d; ++t) + for (int a = 0; a < 3; ++a) + stackFailed += holders_h((k * m.nTri2d + t) * numNodes + a + 3, 1) != + holders_h((next * m.nTri2d + t) * numNodes + a, 1); + } + countFailed += info.dim != 3; + countFailed += info.numVtx != numPlanes * m.nVtx2d; + countFailed += info.numEdge != numLayers * m.nVtx2d; + countFailed += info.numWedge != numWedges; + countFailed += (info.numTri || info.numQuad || info.numTet || info.numHex || + info.numPrism || info.numPyramid); + std::cout << "ids, " << shapeName(ShapeOrder) << ", " + << (periodic ? "periodic" : "open") << ", reverseOrientation " + << reverse << ": id failures " << idFailed << ", range failures " + << rangeFailed << ", stacking failures " << stackFailed + << ", MeshInfo failures " << countFailed << "\n"; + return idFailed + rangeFailed + stackFailed + countFailed; +} + +// Evaluates the field at the nodes and interior points of every wedge of an +// open stack and compares with X. +template +int checkPositions(Omega_h::Mesh &mesh2d, const HostMesh2d &m, Map X, + bool reverse, const char *mapName) { + const auto info = + MeshField::Omegah::getExtrudedWedgeMeshInfo(mesh2d, numLayers); + auto coordField = createCoordField(info); + setExtrudedCoords(coordField.field, m, X); + const auto elem = MeshField::Omegah::getExtrudedWedgeElement( + mesh2d, numLayers, false, reverse); + MeshField::FieldElement fe(info.numWedge, coordField.field, elem.shp, + elem.map); + const auto pts = nodesAndInteriorPoints(); + const int numPts = static_cast(pts.size()); + const auto lc = makeLocalCoords(info.numWedge, pts); + const auto x = MeshField::evaluate(fe, lc, numPts); + const auto x_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x); + int failed = 0; + for (int w = 0; w < info.numWedge; ++w) + for (int p = 0; p < numPts; ++p) { + const Vec3 expected = m.position(X, w, pts[p], reverse); + for (int d = 0; d < 3; ++d) + failed += !isClose(x_h(w * numPts + p, d), expected[d], 1e-12); + } + std::cout << "positions, " << shapeName(ShapeOrder) << ", " << mapName + << ", reverseOrientation " << reverse << ": failures " << failed + << "\n"; + return failed; +} + +// sums w * dV over every point of every wedge +class VolumeIntegrator : public MeshField::Integrator { +public: + VolumeIntegrator(int triOrder, int zetaOrder) + : MeshField::Integrator(triOrder, zetaOrder) {} + void atPoints(Kokkos::View, Kokkos::View w, + Kokkos::View dV) override { + // w holds one weight per point of a single element; dV holds one entry + // per point of every element + const int numPtsPerElem = w.extent(0); + Real sum = 0; + Kokkos::parallel_reduce( + "volume", dV.extent(0), + KOKKOS_LAMBDA(const int i, Real &s) { + s += w(i % numPtsPerElem) * dV(i); + }, + sum); + volume = sum; + } + Real volume = 0; +}; + +// dV is signed, so X must give positively oriented wedges +template +int checkVolume(Omega_h::Mesh &mesh2d, const HostMesh2d &m, Map X, + Real expected, const char *mapName) { + const auto info = + MeshField::Omegah::getExtrudedWedgeMeshInfo(mesh2d, numLayers); + auto coordField = createCoordField(info); + setExtrudedCoords(coordField.field, m, X); + const auto elem = + MeshField::Omegah::getExtrudedWedgeElement(mesh2d, numLayers); + MeshField::FieldElement fe(info.numWedge, coordField.field, elem.shp, + elem.map); + VolumeIntegrator vi(3, 11); + vi.process(fe); + const bool ok = isClose(vi.volume, expected, 1e-12); + std::cout << "volume, " << shapeName(ShapeOrder) << ", " << mapName << ": " + << vi.volume << ", expected " << expected << (ok ? "" : " FAIL") + << "\n"; + return !ok; +} + +// With (R, phi, Z) = (x, s, y), det J has the sign -sign(area) of the 2D +// triangle, and +sign(area) with reverseOrientation. Uses getJacobians, not +// the Integrator, since Debug builds assert dV >= 0 there. +template +int checkOrientation(Omega_h::Mesh &mesh2d, const HostMesh2d &m, bool reverse) { + const auto info = + MeshField::Omegah::getExtrudedWedgeMeshInfo(mesh2d, numLayers); + auto coordField = createCoordField(info); + setExtrudedCoords(coordField.field, m, cylindricalMap); + const auto elem = MeshField::Omegah::getExtrudedWedgeElement( + mesh2d, numLayers, false, reverse); + MeshField::FieldElement fe(info.numWedge, coordField.field, elem.shp, + elem.map); + const int numPts = static_cast(interiorPts.size()); + const auto lc = makeLocalCoords(info.numWedge, interiorPts); + const auto J = MeshField::getJacobians(fe, lc, numPts); + const auto detJ = MeshField::getJacobianDeterminants(fe, J); + const auto detJ_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ); + int failed = 0; + for (int w = 0; w < info.numWedge; ++w) { + const Real area = m.signedArea(w % m.nTri2d); + const Real areaSign = (area > 0) - (area < 0); + const Real expectedSign = reverse ? areaSign : -areaSign; + for (int p = 0; p < numPts; ++p) + failed += !(detJ_h(w * numPts + p) * expectedSign > 0); + } + std::cout << "orientation, " << shapeName(ShapeOrder) + << ", reverseOrientation " << reverse << ": wrong signs " << failed + << " of " << info.numWedge * numPts << "\n"; + return failed; +} + +int runExtrudedWedgeChecks(Omega_h::Mesh &mesh2d) { + const HostMesh2d m(mesh2d); + int failed = 0; + for (bool periodic : {false, true}) + for (bool reverse : {false, true}) { + failed += checkIds<1>(mesh2d, m, periodic, reverse); + failed += checkIds<3>(mesh2d, m, periodic, reverse); + } + for (bool reverse : {false, true}) { + failed += checkPositions<1>(mesh2d, m, taperedMap, reverse, "tapered"); + failed += checkPositions<3>(mesh2d, m, curvedMap, reverse, "curved"); + } + std::cout.precision(16); + failed += checkVolume<1>(mesh2d, m, straightMap, A * H, "straight"); + failed += + checkVolume<1>(mesh2d, m, taperedMap, A * (H + 0.05 * H * H), "tapered"); + failed += checkVolume<3>( + mesh2d, m, curvedMap, + A * (H + 0.05 * H * H + 0.05 * H * H * H + 0.00375 * H * H * H * H), + "sheared and curved"); + for (bool reverse : {false, true}) { + failed += checkOrientation<1>(mesh2d, m, reverse); + failed += checkOrientation<3>(mesh2d, m, reverse); + } + return failed; +} + +} // namespace + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + int numFailed = 0; + { + auto lib = Omega_h::Library(&argc, &argv); + auto mesh2d = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, + 0.0, 2, 2, 0); + numFailed = runExtrudedWedgeChecks(mesh2d); + } + Kokkos::finalize(); + return numFailed ? 1 : 0; +} From d69ce50197b01a74aa0f6b1ca5ac1c62df5c08ff Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Wed, 28 Oct 2026 16:25:16 -0400 Subject: [PATCH 45/45] Clarify extruded wedge numbering --- src/MeshField.hpp | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 9f05e1b..252acd1 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -288,13 +288,17 @@ template auto getTetrahedronElement(Omega_h::Mesh &mesh) { * @brief Maps the nodes of LinearWedgeShape wedges, extruded from a 2D * triangle mesh, to vertices * @details - * The 2D mesh is stacked into numLayers layers and every plane has the same - * triangulation. Vertex v on plane p is vertex p * nVtx2d + v. Wedge w is in - * layer k = w / nTri2d and comes from triangle t = w % nTri2d. Layer k runs - * from plane k to plane k + 1, or to plane (k + 1) % numLayers when periodic. + * The 2D mesh is stacked into numLayers layers. Every plane has the vertex + * numbering and triangle connectivity of the 2D mesh; only the coordinates, + * which the caller sets, differ, so the triangles may change shape from plane + * to plane. Vertex v on plane p is vertex p * nVtx2d + v. Wedge w is in layer + * k = w / nTri2d and comes from triangle t = w % nTri2d. Layer k runs from + * plane k to plane k + 1, or to plane (k + 1) % numLayers when periodic. * Wedge corner a is local vertex a of triangle t, or local vertex {0, 2, 1}[a] * when reverseOrientation is true. Node a is corner a on plane k and node - * a + 3 is corner a on the top plane of the layer. + * a + 3 is corner a on the top plane of the layer. These vertex ids, and the + * vertical edge ids of ExtrudedLinearCubicWedgeToField, number the extruded + * mesh; they are not entity ids of an Omega_h mesh. * * With coordinates (R, phi, Z), phi increasing with the layer and triangles * counterclockwise in (R, Z), det J < 0. Setting reverseOrientation = true