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/31] 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/31] 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/31] 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/31] 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/31] 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/31] 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/31] 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/31] 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/31] change to C-type array --- test/testWedgeElement.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/testWedgeElement.cpp b/test/testWedgeElement.cpp index ac89d46..405f3a2 100644 --- a/test/testWedgeElement.cpp +++ b/test/testWedgeElement.cpp @@ -68,7 +68,7 @@ constexpr int centroidPt = 6; // x(xi) = sum_a L_a(xi0, xi1) [ (1-zeta)/2 B_a + (1+zeta)/2 T_a ] with B_a the // bottom (nodes 0-2) and T_a the top (nodes 3-5) vertex coordinates -Vec3 expectedX(const std::array& nodes, Vec3 xi) { +Vec3 expectedX(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; Vec3 x = Omega_h::zero_vector<3>(); for (int a = 0; a < 3; ++a) { @@ -84,7 +84,7 @@ Vec3 expectedX(const std::array& nodes, Vec3 xi) { // dx/dxi0 = (1-zeta)/2 (B1 - B0) + (1+zeta)/2 (T1 - T0) // dx/dxi1 = (1-zeta)/2 (B2 - B0) + (1+zeta)/2 (T2 - T0) // dx/dzeta = 1/2 sum_a L_a (T_a - B_a) -std::array expectedJ(const std::array& nodes, Vec3 xi) { +std::array expectedJ(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Real down = (1 - xi[2]) / 2; const Real up = (1 + xi[2]) / 2; const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; @@ -170,7 +170,7 @@ int main(int argc, char **argv) { for (int i = 0; i < 3; ++i) for (int j = 0; j < 3; ++j) jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); - Real detExpected = Omega_h::inner_product(jExpected[0] * Omega_h:: cross(jExpected[1], jExpected[2])); + Real detExpected = jExpected[0] * Omega_h:: cross(jExpected[1], jExpected[2]); if (w == 0) detExpected = 2.0 * 3.0 * 4.0 / 2.0; // right prism: a*b*h/2 = 12 detFailed += !isClose(detJ_h(row), detExpected); From 54d471a0c64464a09a1b73fe3f3220df16dfc185 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 15 Oct 2026 14:12:05 -0400 Subject: [PATCH 10/31] rename wedge tests to match class names and split checks into run functions --- ...Element.cpp => testLinearWedgeElement.cpp} | 163 +++++++++--------- ...edgeShape.cpp => testLinearWedgeShape.cpp} | 53 +++--- 2 files changed, 114 insertions(+), 102 deletions(-) rename test/{testWedgeElement.cpp => testLinearWedgeElement.cpp} (54%) rename test/{testWedgeShape.cpp => testLinearWedgeShape.cpp} (80%) diff --git a/test/testWedgeElement.cpp b/test/testLinearWedgeElement.cpp similarity index 54% rename from test/testWedgeElement.cpp rename to test/testLinearWedgeElement.cpp index 405f3a2..0d46d3f 100644 --- a/test/testWedgeElement.cpp +++ b/test/testLinearWedgeElement.cpp @@ -3,9 +3,10 @@ #include "MeshField_Field.hpp" #include "MeshField_For.hpp" #include "MeshField_ShapeField.hpp" -#include #include +#include #include +#include #include // Checks MeshField::LinearWedgeShape inside a MeshField::FieldElement: field @@ -18,15 +19,13 @@ using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; using Vec3 = Omega_h::Vector<3>; - constexpr int numWedges = 2; constexpr int numNodes = MeshField::LinearWedgeShape::numNodes; -bool isClose(Real a, Real b) { - return Omega_h::are_close(a, b, 1e-12, 1e-12); +bool isClose(Real a, Real b) { + return Omega_h::are_close(a, b, 1e-12, 1e-12); } - // Each wedge has its own six vertices: node n of wedge w is vertex 6*w + n. // Nodes 0-2 are the zeta = -1 triangle and nodes 3-5 the zeta = +1 triangle, // matching the LinearWedgeShape node ordering. @@ -63,7 +62,7 @@ const Vec3 localPts[] = { {1. / 3, 1. / 3, 0}, {0.2, 0.3, -0.4}, {0.6, 0.1, 0.7}, {0.05, 0.9, 0.95}}; // clang-format on -constexpr int numPts = sizeof(localPts) / sizeof(localPts[0]); +constexpr int numPts = static_cast(sizeof(localPts) / sizeof(localPts[0])); constexpr int centroidPt = 6; // x(xi) = sum_a L_a(xi0, xi1) [ (1-zeta)/2 B_a + (1+zeta)/2 T_a ] with B_a the @@ -72,11 +71,10 @@ Vec3 expectedX(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Real L[3] = {1 - xi[0] - xi[1], xi[0], xi[1]}; Vec3 x = Omega_h::zero_vector<3>(); for (int a = 0; a < 3; ++a) { - x += L[a] * (1 - xi[2]) / 2 * nodes[a]; - x += L[a] * (1 + xi[2]) / 2 * nodes[a + 3]; + x += L[a] * (1 - xi[2]) / 2 * nodes[a]; + x += L[a] * (1 + xi[2]) / 2 * nodes[a + 3]; } return x; - } // rows of the MeshField Jacobian, J(i,j) = dx_j/dxi_i, from differentiating @@ -91,14 +89,13 @@ std::array expectedJ(const Vec3 (&nodes)[numNodes], Vec3 xi) { const Vec3 dXi0 = down * (nodes[1] - nodes[0]) + up * (nodes[4] - nodes[3]); const Vec3 dXi1 = down * (nodes[2] - nodes[0]) + up * (nodes[5] - nodes[3]); Vec3 dZeta = Omega_h::zero_vector<3>(); - for (int a = 0; a < 3; ++a){ + for (int a = 0; a < 3; ++a) dZeta += L[a] / 2 * (nodes[a + 3] - nodes[a]); - } return {dXi0, dXi1, dZeta}; } - -template +// Writes the wedge vertex coordinates into the coordinate field. +template void setVtxCoords(CoordField field) { const int numVtx = numWedges * numNodes; Kokkos::View coords_h("coords_h", numVtx); @@ -115,74 +112,82 @@ void setVtxCoords(CoordField field) { "setCoordField"); } +// Compares the evaluated positions, Jacobians and determinants of wedge w with +// the hand-differentiated map at every local point. Returns the failure count. +template +int checkWedge(int w, const XView &x_h, const JView &J_h, + const DetView &detJ_h) { + const auto &nodes = wedgeCoords[w]; + int evalFailed = 0, jacFailed = 0, detFailed = 0; + for (int p = 0; p < numPts; ++p) { + const int row = w * numPts + p; + // at a node, x is that node; at the centroid, x is the mean of the six + // vertices; everywhere else x follows the map + Vec3 xExpected = expectedX(nodes, localPts[p]); + if (p < numNodes) + xExpected = nodes[p]; + if (p == centroidPt) { + xExpected = Omega_h::zero_vector<3>(); + for (int n = 0; n < numNodes; ++n) + xExpected += nodes[n] / Real(numNodes); + } + for (int d = 0; d < 3; ++d) + evalFailed += !isClose(x_h(row, d), xExpected[d]); + + const auto jExpected = expectedJ(nodes, localPts[p]); + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); + const Real detExpected = + jExpected[0] * Omega_h::cross(jExpected[1], jExpected[2]); + detFailed += !isClose(detJ_h(row), detExpected); + } + std::cout << "LinearWedgeShape FieldElement, " << wedgeNames[w] + << ": evaluate failures " << evalFailed << ", Jacobian failures " + << jacFailed << ", det J failures " << detFailed << "\n"; + return evalFailed + jacFailed + detFailed; +} + +// Builds the coordinate field and the FieldElement, runs the evaluation, +// Jacobian and determinant, and checks every wedge. Must be called between +// Kokkos::initialize and Kokkos::finalize. +int runLinearWedgeElementChecks() { + const MeshField::MeshInfo meshInfo{ + .numVtx = numWedges * numNodes, .numWedge = numWedges, .dim = 3}; + auto coordField = MeshField::CreateCoordinateField< + ExecutionSpace, MeshField::KokkosController, 3>(meshInfo); + setVtxCoords(coordField.field); + + MeshField::FieldElement fe(numWedges, coordField.field, + MeshField::LinearWedgeShape(), + WedgeToVertexField()); + + // the same local points in every wedge, stored element-major + Kokkos::View lc("localCoords", numWedges * numPts, 3); + auto lc_h = Kokkos::create_mirror_view(lc); + for (int w = 0; w < numWedges; ++w) + for (int p = 0; p < numPts; ++p) + for (int d = 0; d < 3; ++d) + lc_h(w * numPts + p, d) = localPts[p][d]; + Kokkos::deep_copy(lc, lc_h); + + const auto x = MeshField::evaluate(fe, lc, numPts); + const auto J = MeshField::getJacobians(fe, lc, numPts); + const auto detJ = MeshField::getJacobianDeterminants(fe, J); + const auto x_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x); + const auto J_h = Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), J); + const auto detJ_h = + Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ); + + int numFailed = 0; + for (int w = 0; w < numWedges; ++w) + numFailed += checkWedge(w, x_h, J_h, detJ_h); + return numFailed; +} int main(int argc, char **argv) { Kokkos::initialize(argc, argv); - int numFailed = 0; - { - const MeshField::MeshInfo meshInfo{.numVtx = numWedges * numNodes, .numWedge = numWedges, .dim = 3}; - auto coordField = MeshField::CreateCoordinateField< - ExecutionSpace, MeshField::KokkosController, 3>(meshInfo); - setVtxCoords(coordField.field); - - MeshField::FieldElement fe(numWedges, coordField.field, - MeshField::LinearWedgeShape(), - WedgeToVertexField()); - - // the same local points in every wedge, stored element-major - Kokkos::View lc("localCoords", numWedges * numPts, 3); - auto lc_h = Kokkos::create_mirror_view(lc); - for (int w = 0; w < numWedges; ++w) - for (int p = 0; p < numPts; ++p) - for (int d = 0; d < 3; ++d) - lc_h(w * numPts + p, d) = localPts[p][d]; - Kokkos::deep_copy(lc, lc_h); - - const auto x = MeshField::evaluate(fe, lc, numPts); - const auto J = MeshField::getJacobians(fe, lc, numPts); - const auto detJ = MeshField::getJacobianDeterminants(fe, J); - const auto x_h = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), x); - const auto J_h = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), J); - const auto detJ_h = - Kokkos::create_mirror_view_and_copy(Kokkos::HostSpace(), detJ); - - for (int w = 0; w < numWedges; ++w) { - const auto &nodes = wedgeCoords[w]; - int evalFailed = 0, jacFailed = 0, detFailed = 0; - for (int p = 0; p < numPts; ++p) { - const int row = w * numPts + p; - // evaluate: x at a node is that node, x at the centroid is the mean - // of the six vertices, and x matches the map x(xi) everywhere - Vec3 xExpected = expectedX(nodes, localPts[p]); - if (p < numNodes) - xExpected = nodes[p]; - if (p == centroidPt) { - xExpected = Omega_h::zero_vector<3>(); - for (int n = 0; n < numNodes; ++n) - xExpected += nodes[n]/numNodes ; - } - for (int d = 0; d < 3; ++d) - evalFailed += !isClose(x_h(row, d), xExpected[d]); - // Jacobian rows and determinant - const auto jExpected = expectedJ(nodes, localPts[p]); - for (int i = 0; i < 3; ++i) - for (int j = 0; j < 3; ++j) - jacFailed += !isClose(J_h(row, i, j), jExpected[i][j]); - Real detExpected = jExpected[0] * Omega_h:: cross(jExpected[1], jExpected[2]); - if (w == 0) - detExpected = 2.0 * 3.0 * 4.0 / 2.0; // right prism: a*b*h/2 = 12 - detFailed += !isClose(detJ_h(row), detExpected); - } - std::cout << "LinearWedgeShape FieldElement, " << wedgeNames[w] - << ": evaluate failures " << evalFailed - - << ", Jacobian failures " << jacFailed << ", det J failures " - << detFailed << "\n"; - numFailed += evalFailed + jacFailed + detFailed; - } - } + const int numFailed = runLinearWedgeElementChecks(); Kokkos::finalize(); return numFailed ? 1 : 0; -} +} \ No newline at end of file diff --git a/test/testWedgeShape.cpp b/test/testLinearWedgeShape.cpp similarity index 80% rename from test/testWedgeShape.cpp rename to test/testLinearWedgeShape.cpp index b5a8aa6..267c966 100644 --- a/test/testWedgeShape.cpp +++ b/test/testLinearWedgeShape.cpp @@ -1,18 +1,17 @@ -#include +#include #include #include #include // Checks of MeshField::LinearWedgeShape. The checks run inside Kokkos kernels -// so the shape functions are exercised in the default execution space -// and failures are counted explicitly so the test also fails in -// builds where assert() is disabled. +// so the shape functions are exercised in the default execution space, and +// failures are counted explicitly so the test also fails in builds where +// assert() is disabled. using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MeshField::Real; using Shape = MeshField::LinearWedgeShape; - KOKKOS_INLINE_FUNCTION bool isClose(Real a, Real b, Real tol) { return Kokkos::fabs(a - b) <= tol; } @@ -46,10 +45,8 @@ int checkKroneckerDelta() { return numFailed; } -// At interior points: the values sum to one, the gradients sum to zero, the -// basis reproduces f and its gradient exactly, and the analytic gradients -// match central finite differences of the values. -int checkInteriorPoints() { +// Interior points, spanning the zeta range. +Kokkos::View makeInteriorPoints() { // clang-format off const Real pts[][3] = {{1. / 3, 1. / 3, 0.0}, {0.2, 0.3, -0.4}, @@ -57,13 +54,22 @@ int checkInteriorPoints() { {0.05, 0.9, 0.95}, {0.45, 0.45, -0.95}}; // clang-format on - const int numPts = sizeof(pts) / sizeof(pts[0]); + const int numPts = static_cast(sizeof(pts) / sizeof(pts[0])); Kokkos::View pts_d("wedgePts", numPts); auto pts_h = Kokkos::create_mirror_view(pts_d); for (int i = 0; i < numPts; ++i) for (int d = 0; d < 3; ++d) pts_h(i, d) = pts[i][d]; Kokkos::deep_copy(pts_d, pts_h); + return pts_d; +} + +// At interior points: the values sum to one, the gradients sum to zero, the +// basis reproduces f and its gradient exactly, and the analytic gradients +// match central finite differences of the values. +int checkInteriorPoints() { + const auto pts_d = makeInteriorPoints(); + const int numPts = static_cast(pts_d.extent(0)); int numFailed = 0; Kokkos::parallel_reduce( @@ -97,8 +103,7 @@ int checkInteriorPoints() { failed += !isClose(gInterp[d], gExact[d], tol); } - // gradient check with finite difference evaluation - + // gradients against central finite differences of the values const Real h = 1e-6; for (int d = 0; d < 3; ++d) { auto xp = xi; @@ -115,19 +120,21 @@ int checkInteriorPoints() { return numFailed; } +// Runs both checks and prints their failure counts. Must be called between +// Kokkos::initialize and Kokkos::finalize. +int runLinearWedgeShapeChecks() { + const int kronecker = checkKroneckerDelta(); + std::cout << "LinearWedgeShape Kronecker delta failures: " << kronecker + << "\n"; + const int interior = checkInteriorPoints(); + std::cout << "LinearWedgeShape interior point failures: " << interior + << "\n"; + return kronecker + interior; +} int main(int argc, char **argv) { Kokkos::initialize(argc, argv); - int numFailed = 0; - { - const int kronecker = checkKroneckerDelta(); - std::cout << "LinearWedgeShape Kronecker delta failures: " << kronecker - << "\n"; - const int interior = checkInteriorPoints(); - std::cout << "LinearWedgeShape interior point failures: " << interior - << "\n"; - numFailed = kronecker + interior; - } + const int numFailed = runLinearWedgeShapeChecks(); Kokkos::finalize(); return numFailed ? 1 : 0; -} +} \ No newline at end of file From a198061e36a77877dd6734d2fc946a1d40a96ab0 Mon Sep 17 00:00:00 2001 From: Abhiyan Paudel <83005772+abhiyanpaudel@users.noreply.github.com> Date: Thu, 15 Oct 2026 14:18:15 -0400 Subject: [PATCH 11/31] 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 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 12/31] 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 13/31] 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 14/31] 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 15/31] 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 16/31] 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 17/31] 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 18/31] 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 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 19/31] 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 20/31] 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 21/31] 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 22/31] 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 23/31] 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 24/31] 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 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 25/31] 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 26/31] 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 27/31] 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 28/31] 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 29/31] 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 30/31] 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 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 31/31] 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() {