From 19b9334cf0fb9b572bcf072a50fd4a7ef86a427e Mon Sep 17 00:00:00 2001 From: Sichao25 Date: Fri, 12 Jun 2026 02:30:00 -0400 Subject: [PATCH 1/5] create reduced quintic element backend --- CMakeLists.txt | 5 + src/MeshField.hpp | 132 ++++++ src/MeshField_Element.hpp | 84 +++- src/MeshField_ReducedQuintic.hpp | 446 ++++++++++++++++++ src/MeshField_Shape.hpp | 259 +++++++++++ test/testOmegahTriReducedQuintic.cpp | 305 ++++++++++++ test/testReducedQuintic.cpp | 663 +++++++++++++++++++++++++++ 7 files changed, 1885 insertions(+), 9 deletions(-) create mode 100644 src/MeshField_ReducedQuintic.hpp create mode 100644 test/testOmegahTriReducedQuintic.cpp create mode 100644 test/testReducedQuintic.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 4642c49..f4402af 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -36,6 +36,7 @@ set(MESHFIELD_HEADERS src/MeshField_Field.hpp src/MeshField.hpp src/MeshField_SPR_ErrorEstimator.hpp + src/MeshField_ReducedQuintic.hpp "${CMAKE_CURRENT_BINARY_DIR}/MeshField_Config.hpp" ) if(MeshFields_USE_Cabana) @@ -198,6 +199,8 @@ meshfields_add_exe(ExceptionTest test/testExceptions.cpp) meshfields_add_exe(PointMapping test/testPointMapping.cpp) meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) +meshfields_add_exe(ReducedQuinticTest test/testReducedQuintic.cpp) +meshfields_add_exe(OmegahReducedQuinticTest test/testOmegahTriReducedQuintic.cpp) if(MeshFields_USE_Cabana) meshfields_add_exe(ControllerPerformance test/testControllerPerformance.cpp) @@ -216,6 +219,8 @@ test_func(CountIntegrator ./CountIntegrator) smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(PointMapping ./PointMapping) test_func(OmegahTetTest ./OmegahTetTest) +test_func(ReducedQuinticTest ./ReducedQuinticTest) +test_func(OmegahReducedQuinticTest ./OmegahReducedQuinticTest) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 2424710..5def3d9 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -6,9 +6,11 @@ #include "MeshField_Fail.hpp" #include "MeshField_For.hpp" #include "MeshField_ShapeField.hpp" +#include "MeshField_ReducedQuintic.hpp" #include "Omega_h_file.hpp" //move #include "Omega_h_mesh.hpp" //move #include "Omega_h_simplex.hpp" //move +#include namespace { @@ -238,6 +240,65 @@ struct QuadraticTetrahedronToField { } }; +struct ReducedQuinticTriangleToField { + Omega_h::LOs triVerts; + + ReducedQuinticTriangleToField(Omega_h::Mesh& mesh) + : triVerts(mesh.ask_elem_verts()) { + + if (mesh.dim() != 2 || mesh.family() != OMEGA_H_SIMPLEX) { + MeshField::fail( + "The mesh passed to %s must be 2D and simplex (triangles)\n", + __func__); + } + } + + static constexpr KOKKOS_FUNCTION + Kokkos::Array + getTopology() { + return {MeshField::Triangle}; + } + + KOKKOS_FUNCTION + MeshField::ElementToDofHolderMap + operator()(MeshField::LO triNodeIdx, + MeshField::LO triCompIdx, + MeshField::LO tri, + MeshField::Mesh_Topology topo) const { + + assert(topo == MeshField::Triangle); + + // shape function index -> vertex + const MeshField::LO localVtxIdx = triNodeIdx / 6; + + // dof component within vertex + const MeshField::LO vertexDof = triNodeIdx % 6; + + const auto triDim = 2; + const auto vtxDim = 0; + const auto ignored = -1; + + const auto canonicalVtxIdx = + (Omega_h::simplex_down_template( + triDim, vtxDim, localVtxIdx, ignored) + + 2) % + 3; + + const auto triToVtxDegree = + Omega_h::simplex_degree(triDim, vtxDim); + + const MeshField::LO vtx = + triVerts[(tri * triToVtxDegree) + canonicalVtxIdx]; + + return { + vertexDof, // node within vertex dof holder + triCompIdx, // field component + vtx, // vertex entity + MeshField::Vertex // topology + }; + } +}; + template auto getTriangleElement(Omega_h::Mesh &mesh) { static_assert(ShapeOrder == 1 || ShapeOrder == 2); if constexpr (ShapeOrder == 1) { @@ -256,6 +317,50 @@ template auto getTriangleElement(Omega_h::Mesh &mesh) { QuadraticTriangleToField(mesh)}; } } + +// ReducedQuintic triangle element with precomputed coefficients +inline auto getReducedQuinticTriangleElement(Omega_h::Mesh &mesh) { + if (mesh.dim() != 2 || mesh.family() != OMEGA_H_SIMPLEX) { + MeshField::fail("getReducedQuinticTriangleElement requires 2D simplex mesh\n"); + } + + const auto numTri = mesh.nfaces(); + const auto coords_d = mesh.coords(); + const auto triVerts_d = mesh.ask_elem_verts(); + + Omega_h::HostRead triVerts(triVerts_d); + Omega_h::HostRead coords(coords_d); + + // Extract triangle vertex coordinates + std::vector triCoords(numTri * 6); // 3 vertices * 2 coords per triangle + + for (Omega_h::LO tri = 0; tri < numTri; tri++) { + for (int vi = 0; vi < 3; vi++) { + const auto triDim = 2; + const auto vtxDim = 0; + const auto ignored = -1; + const auto localVtxIdx = (Omega_h::simplex_down_template(triDim, vtxDim, vi, ignored) + 2) % 3; + const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim); + const Omega_h::LO vtx = triVerts[tri * triToVtxDegree + localVtxIdx]; + + triCoords[tri * 6 + vi * 2 + 0] = coords[vtx * 2 + 0]; // x + triCoords[tri * 6 + vi * 2 + 1] = coords[vtx * 2 + 1]; // y + } + } + + // Precompute coefficients for all triangles + auto elemCoeffs = MeshField::precomputeReducedQuinticCoefficients(numTri, triCoords.data()); + + struct result { + MeshField::ReducedQuinticTriangleShape shp; + ReducedQuinticTriangleToField map; + Kokkos::View coeffs; + }; + + return result{MeshField::ReducedQuinticTriangleShape(), + ReducedQuinticTriangleToField(mesh), + elemCoeffs}; +} template auto getTetrahedronElement(Omega_h::Mesh &mesh) { static_assert(ShapeOrder == 1 || ShapeOrder == 2); if constexpr (ShapeOrder == 1) { @@ -360,6 +465,33 @@ class OmegahMeshField { return eval; } + // evaluate a ReducedQuintic field at the specified local coordinates for each triangle + template + auto triangleReducedQuinticEval(const ViewType &localCoords, + Kokkos::View offsets, + const ShapeField &field) const { + const auto MeshDim = 2; + if (mesh.dim() != MeshDim) { + MeshField::fail("input mesh must be 2d\n"); + } + + const auto [shp, map, coeffs] = Omegah::getReducedQuinticTriangleElement(mesh); + + MeshField::FieldElement f( + meshInfo.numTri, field, shp, map, coeffs); + auto eval = MeshField::evaluate(f, localCoords, offsets); + return eval; + } + + // evaluate a ReducedQuintic field at the specified local coordinate for each triangle + template + auto triangleReducedQuinticEval(const ViewType &localCoords, + size_t NumPtsPerElem, + const ShapeField &field) const { + auto offsets = createOffsets(meshInfo.numTri, NumPtsPerElem); + return triangleReducedQuinticEval(localCoords, offsets, field); + } + template auto tetrahedronLocalPointEval(const ViewType &localCoords, size_t NumPtsPerElem, diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index dfe302f..24678a1 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -5,6 +5,7 @@ #include #include #include +#include #include // getLastValue #include #include // has_static_size helper @@ -130,13 +131,24 @@ struct FieldElement { const FieldAccessor field; ShapeType shapeFn; ElementDofHolderAccessor elm2dof; + Kokkos::View elemCoeffs; // Optional: per-element coefficients for advanced shapes static const size_t MeshEntDim = ShapeType::meshEntDim; + + // Standard constructor (for shapes without per-element coefficients) FieldElement(size_t numMeshEntsIn, const FieldAccessor &fieldIn, const ShapeType shapeFnIn, const ElementDofHolderAccessor elm2dofIn) : numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn), - elm2dof(elm2dofIn) {} + elm2dof(elm2dofIn), elemCoeffs() {} + + // Constructor with per-element coefficients (for reduced quintic, etc.) + FieldElement(size_t numMeshEntsIn, const FieldAccessor &fieldIn, + const ShapeType shapeFnIn, + const ElementDofHolderAccessor elm2dofIn, + Kokkos::View elemCoeffsIn) + : numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn), + elm2dof(elm2dofIn), elemCoeffs(elemCoeffsIn) {} /* general template for baseType which simply sets type */ template struct baseType { @@ -173,16 +185,70 @@ struct FieldElement { assert(ent >= 0); assert(static_cast(ent) < numMeshEnts); ValArray c; - const auto shapeValues = shapeFn.getValues(localCoord); + + // Get shape function values + // For shapes that require per-element coefficients (e.g., reduced quintic), pass the coefficients to the shape function + auto shapeValues = [&]() { + if constexpr (std::is_same_v, + ReducedQuinticTriangleShape>) { + assert(elemCoeffs.data() != nullptr); + auto coeffSlice = Kokkos::subview(elemCoeffs, ent, Kokkos::ALL()); + return shapeFn.getValues(localCoord, coeffSlice); + } else { + return shapeFn.getValues(localCoord); + } + }(); + for (size_t ci = 0; ci < NumComponents; ++ci) c[ci] = 0; - for (auto topo : elm2dof.getTopology()) { // element topology - for (size_t ni = 0; ni < shapeFn.numNodes; ++ni) { - for (size_t ci = 0; ci < NumComponents; ++ci) { - auto map = elm2dof(ni, ci, ent, topo); - const auto fval = - field(map.entity, map.node, map.component, map.topo); - c[ci] += fval * shapeValues[ni]; + + // For ReducedQuintic, we need to transform DOFs from physical to local coordinates + if constexpr (std::is_same_v, + ReducedQuinticTriangleShape>) { + assert(elemCoeffs.data() != nullptr); + auto coeffSlice = Kokkos::subview(elemCoeffs, ent, Kokkos::ALL()); + + // Extract rotation parameters from coefficients + // elemCoeffs layout: [order[0], order[1], order[2], a, b, c, sin_theta, cos_theta, ...] + const Real sin_theta = coeffSlice(6); + const Real cos_theta = coeffSlice(7); + + for (auto topo : elm2dof.getTopology()) { + // ReducedQuintic has 18 nodes: 3 vertices × 6 DOFs per vertex + const size_t numVertices = 3; + const size_t dofsPerVertex = 6; + + for (size_t vi = 0; vi < numVertices; ++vi) { + // Gather all 6 DOFs for this vertex in physical coordinates + Real physicalDofs[6]; + for (size_t di = 0; di < dofsPerVertex; ++di) { + const size_t ni = vi * dofsPerVertex + di; + auto map = elm2dof(ni, 0, ent, topo); // ci=0 since ReducedQuintic is scalar + physicalDofs[di] = field(map.entity, map.node, map.component, map.topo); + } + + // Transform DOFs to local coordinates and accumulate + for (size_t di = 0; di < dofsPerVertex; ++di) { + const size_t ni = vi * dofsPerVertex + di; + const Real localDof = transformDofPhysicalToLocal( + di, sin_theta, cos_theta, physicalDofs); + + for (size_t ci = 0; ci < NumComponents; ++ci) { + c[ci] += localDof * shapeValues[ni]; + } + } + } + } + } else { + // Standard DOF gathering for other shape functions + for (auto topo : elm2dof.getTopology()) { + for (size_t ni = 0; ni < shapeFn.numNodes; ++ni) { + for (size_t ci = 0; ci < NumComponents; ++ci) { + auto map = elm2dof(ni, ci, ent, topo); + const auto fval = + field(map.entity, map.node, map.component, map.topo); + c[ci] += fval * shapeValues[ni]; + } } } } diff --git a/src/MeshField_ReducedQuintic.hpp b/src/MeshField_ReducedQuintic.hpp new file mode 100644 index 0000000..5999126 --- /dev/null +++ b/src/MeshField_ReducedQuintic.hpp @@ -0,0 +1,446 @@ +#ifndef MESHFIELD_REDUCEDQUINTIC_HPP +#define MESHFIELD_REDUCEDQUINTIC_HPP + +#include +#include +#include +#include +#include +#include + +namespace MeshField { + +/** + * @brief Simple LU decomposition with partial pivoting (internal helper) + * + * Solves the linear system A*X = B where A is n×n and B is n×nrhs. + * This replaces the need for LAPACK's dgesv for our 20×20 system. + * Works with flattened arrays internally for the actual computation. + * + * @param n Dimension of the matrix (20 for reduced quintic) + * @param nrhs Number of right-hand sides (18 for reduced quintic) + * @param A Input matrix (destroyed on output), row-major format + * @param lda Leading dimension of A + * @param B Input/output: RHS on input, solution on output, row-major format + * @param ldb Leading dimension of B + * @return 0 on success, positive value if singular + */ +inline int solveLU_internal(int n, int nrhs, Real* A, int lda, Real* B, int ldb) { + std::vector pivot(n); + const Real eps = 1e-14; + + // LU decomposition with partial pivoting + for (int k = 0; k < n; k++) { + // Find pivot + int maxRow = k; + Real maxVal = std::abs(A[k * lda + k]); + for (int i = k + 1; i < n; i++) { + Real val = std::abs(A[i * lda + k]); + if (val > maxVal) { + maxVal = val; + maxRow = i; + } + } + + if (maxVal < eps) { + return k + 1; // Singular matrix + } + + pivot[k] = maxRow; + + // Swap rows in A if needed + if (maxRow != k) { + for (int j = 0; j < n; j++) { + std::swap(A[k * lda + j], A[maxRow * lda + j]); + } + } + + // Eliminate column + Real diag = A[k * lda + k]; + for (int i = k + 1; i < n; i++) { + Real factor = A[i * lda + k] / diag; + A[i * lda + k] = factor; // Store multiplier + for (int j = k + 1; j < n; j++) { + A[i * lda + j] -= factor * A[k * lda + j]; + } + } + } + + // Apply pivoting to B + for (int k = 0; k < n; k++) { + if (pivot[k] != k) { + for (int rhs = 0; rhs < nrhs; rhs++) { + std::swap(B[k * ldb + rhs], B[pivot[k] * ldb + rhs]); + } + } + } + + // Forward substitution: solve L*Y = B for each column + for (int i = 0; i < n; i++) { + for (int j = 0; j < i; j++) { + Real factor = A[i * lda + j]; + for (int rhs = 0; rhs < nrhs; rhs++) { + B[i * ldb + rhs] -= factor * B[j * ldb + rhs]; + } + } + } + + // Backward substitution: solve U*X = Y for each column + for (int i = n - 1; i >= 0; i--) { + Real diag = A[i * lda + i]; + for (int j = i + 1; j < n; j++) { + Real factor = A[i * lda + j]; + for (int rhs = 0; rhs < nrhs; rhs++) { + B[i * ldb + rhs] -= factor * B[j * ldb + rhs]; + } + } + for (int rhs = 0; rhs < nrhs; rhs++) { + B[i * ldb + rhs] /= diag; + } + } + + return 0; // Success +} + +/** + * @brief Rotate DOFs for reduced quintic element + * + * @param dofs_p Input DOFs in original orientation (6 DOFs per vertex) + * @param sin_theta_p Sine of rotation angle + * @param cos_theta_p Cosine of rotation angle + */ +void rotateDof(double dofs_p[6], double sin_theta_p, double cos_theta_p) +{ + const double s = sin_theta_p; + const double c = cos_theta_p; + const double ss = s * s; + const double cc = c * c; + const double sc = s * c; + + const double dofs_r[6] = { + dofs_p[0], + c * dofs_p[1] + s * dofs_p[2], + c * dofs_p[2] - s * dofs_p[1], + cc * dofs_p[3] + 2*sc * dofs_p[4] + ss * dofs_p[5], + -sc * dofs_p[3] + (cc - ss) * dofs_p[4] + sc * dofs_p[5], + ss * dofs_p[3] - 2*sc * dofs_p[4] + cc * dofs_p[5] + }; + + for (int i = 0; i < 6; i++) + dofs_p[i] = dofs_r[i]; +} + +/** + * @brief Device-compatible function to transform DOFs from physical to local coordinates + * + * Transforms derivatives from physical (x,y) coordinates to local (ξ,η) coordinates + * using a rotation. This is necessary because ReducedQuintic shape functions are + * defined in a local coordinate system aligned with the triangle's geometry. + * + * @param dof_value The DOF value to transform (one of the 6 components per vertex) + * @param vertex_idx Vertex index (0, 1, or 2) + * @param dof_idx DOF component index (0-5: value, ∂x, ∂y, ∂²x², ∂²xy, ∂²y²) + * @param sin_theta Sine of rotation angle (from elemCoeffs) + * @param cos_theta Cosine of rotation angle (from elemCoeffs) + * @param allDofs Array of all 6 DOFs for this vertex in physical coordinates + * @return The transformed DOF value in local coordinates + */ +template +KOKKOS_INLINE_FUNCTION +Real transformDofPhysicalToLocal( + int dof_idx, + Real sin_theta, + Real cos_theta, + const Real allDofs[6]) +{ + const Real s = sin_theta; + const Real c = cos_theta; + const Real ss = s * s; + const Real cc = c * c; + const Real sc = s * c; + + // Transform based on DOF type + switch (dof_idx) { + case 0: // value - unchanged + return allDofs[0]; + case 1: // ∂/∂x → ∂/∂ξ + return c * allDofs[1] + s * allDofs[2]; + case 2: // ∂/∂y → ∂/∂η + return c * allDofs[2] - s * allDofs[1]; + case 3: // ∂²/∂x² + return cc * allDofs[3] + 2*sc * allDofs[4] + ss * allDofs[5]; + case 4: // ∂²/∂x∂y + return -sc * allDofs[3] + (cc - ss) * allDofs[4] + sc * allDofs[5]; + case 5: // ∂²/∂y² + return ss * allDofs[3] - 2*sc * allDofs[4] + cc * allDofs[5]; + default: + return allDofs[dof_idx]; + } +} + +/** + * @brief Reorder triangle vertices to put longest edge along local xi-axis + * + * Reordering strategy: + * - Find the longest edge + * - Orient it to be the local xi-axis (from v0 to v1) + * - Ensure counter-clockwise ordering + * + * @param coords Triangle vertex coordinates [3][2] + * @param order Output vertex order [3] - order[i] gives original index of reordered vertex i + */ +inline void reorderTriangleVertices( + const Real coords[3][2], + int order[3]) +{ + // Compute edge vectors + Real v1v2[2] = {coords[1][0] - coords[0][0], coords[1][1] - coords[0][1]}; + Real v1v3[2] = {coords[2][0] - coords[0][0], coords[2][1] - coords[0][1]}; + Real v2v3[2] = {coords[2][0] - coords[1][0], coords[2][1] - coords[1][1]}; + + // Check if counter-clockwise + int counterclockwise = 1; + if (v1v2[0]*v1v3[1] - v1v2[1]*v1v3[0] < 0) { + counterclockwise = 0; + } + + // Compute edge length squares + Real v1v2lensq = v1v2[0]*v1v2[0] + v1v2[1]*v1v2[1]; + Real v1v3lensq = v1v3[0]*v1v3[0] + v1v3[1]*v1v3[1]; + Real v2v3lensq = v2v3[0]*v2v3[0] + v2v3[1]*v2v3[1]; + + // Find longest edge and set order + if (v1v2lensq > v1v3lensq && v1v2lensq > v2v3lensq) { + // Edge v0-v1 is longest + if (counterclockwise) { + order[0] = 0; + order[1] = 1; + order[2] = 2; + } else { + order[0] = 1; + order[1] = 0; + order[2] = 2; + } + } else { + if (v1v3lensq > v2v3lensq) { + // Edge v0-v2 is longest + if (counterclockwise) { + order[0] = 2; + order[1] = 0; + order[2] = 1; + } else { + order[0] = 0; + order[1] = 2; + order[2] = 1; + } + } else { + // Edge v1-v2 is longest + if (counterclockwise) { + order[0] = 1; + order[1] = 2; + order[2] = 0; + } else { + order[0] = 2; + order[1] = 1; + order[2] = 0; + } + } + } +} + +/** + * @brief Compute geometric parameters for reduced quintic triangle element + * + * Computes geometric parameters by: + * 1. Reordering vertices to put longest edge along local xi-axis + * 2. Computing local coordinate system parameters + * + * @param coords Original triangle vertex coordinates [3][2] + * @param origin Output: origin of local coordinate system + * @param a Output: distance from origin to first vertex in local system + * @param b Output: distance from origin to second vertex in local system + * @param c Output: perpendicular distance to third vertex + * @param sin_theta Output: sine of rotation angle + * @param cos_theta Output: cosine of rotation angle + * @param order Output: vertex reordering [3] - order[i] gives original index + */ +inline void computeReducedQuinticGeometry( + const Real coords[3][2], + Real origin[2], + Real& a, Real& b, Real& c, + Real& sin_theta, Real& cos_theta, + int order[3]) +{ + // First reorder vertices to put longest edge along xi-axis + reorderTriangleVertices(coords, order); + + // Apply reordering + Real abcCoord[3][2]; + for (int i = 0; i < 3; i++) { + int idx = order[i]; + abcCoord[i][0] = coords[idx][0]; + abcCoord[i][1] = coords[idx][1]; + } + + // Coordinate system (after reordering): + // - Origin at projection of v2 onto v0-v1 edge + // - v0 at (-b, 0) relative to origin + // - v1 at (a, 0) relative to origin + // - v2 at (0, c) relative to origin + + // Edge from v0 to v1 (now the longest edge) + Real ab[2] = {abcCoord[1][0] - abcCoord[0][0], abcCoord[1][1] - abcCoord[0][1]}; + Real ablensq = ab[0]*ab[0] + ab[1]*ab[1]; + Real ablen = std::sqrt(ablensq); + + // Vector from v0 to v2 + Real ac[2] = {abcCoord[2][0] - abcCoord[0][0], abcCoord[2][1] - abcCoord[0][1]}; + + // Project v2 onto v0-v1 edge to find origin + Real buffer = (ab[0]*ac[0] + ab[1]*ac[1]) / ablensq; + Real ac2ab[2] = {ab[0] * buffer, ab[1] * buffer}; + origin[0] = ac2ab[0] + abcCoord[0][0]; + origin[1] = ac2ab[1] + abcCoord[0][1]; + + // Distances from origin to each vertex (in reordered system) + Real vec_a[2] = {origin[0] - abcCoord[0][0], origin[1] - abcCoord[0][1]}; + Real vec_b[2] = {origin[0] - abcCoord[1][0], origin[1] - abcCoord[1][1]}; + Real vec_c[2] = {origin[0] - abcCoord[2][0], origin[1] - abcCoord[2][1]}; + + b = std::sqrt(vec_a[0]*vec_a[0] + vec_a[1]*vec_a[1]); // distance to v0 + a = std::sqrt(vec_b[0]*vec_b[0] + vec_b[1]*vec_b[1]); // distance to v1 + c = std::sqrt(vec_c[0]*vec_c[0] + vec_c[1]*vec_c[1]); // distance to v2 + + // Verify: a + b should equal edge length + assert(std::abs(a + b - ablen) < 1e-10 * ablen); + + // Angle of v0-v1 edge + sin_theta = ab[1] / ablen; + cos_theta = ab[0] / ablen; +} + +/** + * @brief Compute coefficient matrix for reduced quintic element + * + * @param a, b, c Geometric parameters from computeReducedQuinticGeometry + * @param coeffs Output: coefficient matrix [18][20] + */ +inline void computeReducedQuinticCoefficients( + Real a, Real b, Real c, + Real coeffs[18][20]) +{ + // Compute powers + const Real a2 = a*a, a3 = a2*a, a4 = a2*a2, a5 = a2*a3; + const Real b2 = b*b, b3 = b2*b, b4 = b2*b2, b5 = b2*b3; + const Real c2 = c*c, c3 = c2*c, c4 = c2*c2, c5 = c2*c3; + + // Build constraint matrix A (20x20) + // This encodes the boundary conditions at vertices and continuity constraints + Real A[20][20] = { + {1., -b, 0., b2, 0, 0, -b3, 0, 0, 0, b4, 0, 0, 0, 0, -b5, 0, 0, 0, 0}, + {0., 1., 0., -2.*b, 0., 0., 3*b2, 0, 0, 0, -4*b3, 0, 0, 0, 0, 5*b4, 0, 0, 0, 0}, + {0, 0, 1, 0, -b, 0, 0, b2, 0, 0, 0, -b3, 0, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 2, 0, 0, -6*b, 0, 0, 0, 12*b2, 0, 0, 0, 0, -20*b3, 0, 0, 0, 0}, + {0, 0, 0, 0, 1, 0, 0, -2*b, 0, 0, 0, 3*b2, 0, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 0, 2, 0, 0, -2*b, 0, 0, 0, 2*b2, 0, 0, 0, -2*b3, 0, 0, 0}, + {1, a, 0, a2, 0, 0, a3, 0, 0, 0, a4, 0, 0, 0, 0, a5, 0, 0, 0, 0}, + {0, 1, 0, 2*a, 0, 0, 3*a2, 0, 0, 0, 4*a3, 0, 0, 0, 0, 5*a4, 0, 0, 0, 0}, + {0, 0, 1, 0, a, 0, 0, a2, 0, 0, 0, a3, 0, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 2, 0, 0, 6*a, 0, 0, 0, 12*a2, 0, 0, 0, 0, 20*a3, 0, 0, 0, 0}, + {0, 0, 0, 0, 1, 0, 0, 2*a, 0, 0, 0, 3*a2, 0, 0, 0, 0, 0, 0, 0, 0}, + {0, 0, 0, 0, 0, 2, 0, 0, 2*a, 0, 0, 0, 2*a2, 0, 0, 0, 2*a3, 0, 0, 0}, + {1, 0, c, 0, 0, c2, 0, 0, 0, c3, 0, 0, 0, 0, c4, 0, 0, 0, 0, c5}, + {0, 1, 0, 0, c, 0, 0, 0, c2, 0, 0, 0, 0, c3, 0, 0, 0, 0, c4, 0}, + {0, 0, 1, 0, 0, 2*c, 0, 0, 0, 3*c2, 0, 0, 0, 0, 4*c3, 0, 0, 0, 0, 5*c4}, + {0, 0, 0, 2, 0, 0, 0, 2*c, 0, 0, 0, 0, 2*c2, 0, 0, 0, 0, 2*c3, 0, 0}, + {0, 0, 0, 0, 1, 0, 0, 0, 2*c, 0, 0, 0, 0, 3*c2, 0, 0, 0, 0, 4*c3, 0}, + {0, 0, 0, 0, 0, 2, 0, 0, 0, 6*c, 0, 0, 0, 0, 12*c2, 0, 0, 0, 0, 20*c3}, + {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5*a4*c, 3*a2*c3-2*a4*c, -2*a*c4+3*a3*c2, c5-4*a2*c3, 5*a*c4}, + {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5*b4*c, 3*b2*c3-2*b4*c, 2*b*c4-3*b3*c2, c5-4*b2*c3, -5*b*c4} + }; + + // Initialize RHS as 20×18 identity matrix (first 18 rows are identity) + Real B[20][18]; + for (int i = 0; i < 20; i++) { + for (int j = 0; j < 18; j++) { + B[i][j] = (i == j) ? 1.0 : 0.0; + } + } + + // Solve A*X = B using our LU solver (row-major) + int info = solveLU_internal(20, 18, &A[0][0], 20, &B[0][0], 18); + + if (info != 0) { + fail("LU solver failed with info = %d in reduced quintic coefficient computation\n", info); + } + + // Copy solution to coeffs[18][20] + for (int j = 0; j < 18; j++) { + for (int i = 0; i < 20; i++) { + coeffs[j][i] = B[i][j]; + } + } +} + +/** + * @brief Precompute all reduced quintic coefficients for a mesh + * + * @param numTriangles Number of triangles in the mesh + * @param triangleCoords Triangle vertex coordinates, shape [numTriangles][3][2] + * @return Kokkos::View with coefficients, shape [numTriangles][6 + 18*20] + * Each row: [order[0], order[1], order[2], a, b, c, coeff_0_0, coeff_0_1, ..., coeff_17_19] + */ +inline Kokkos::View precomputeReducedQuinticCoefficients( + int numTriangles, + const Real* triangleCoords) +{ + // Allocate device view with extended storage for sin_theta and cos_theta + // Layout: [order[0], order[1], order[2], a, b, c, sin_theta, cos_theta, coeff_0_0, ..., coeff_17_19] + Kokkos::View coeffs_d("coefficients_device", numTriangles, 8 + 18 * 20); + + // Create a mirror view on host with matching layout + auto coeffs_h = Kokkos::create_mirror_view(coeffs_d); + + // Compute for each triangle + for (int tri = 0; tri < numTriangles; tri++) { + Real coords[3][2]; + for (int v = 0; v < 3; v++) { + coords[v][0] = triangleCoords[tri * 6 + v * 2 + 0]; + coords[v][1] = triangleCoords[tri * 6 + v * 2 + 1]; + } + + // Compute geometric parameters + Real origin[2], a, b, c, sin_theta, cos_theta; + int order[3]; + computeReducedQuinticGeometry(coords, origin, a, b, c, sin_theta, cos_theta, order); + + // Store vertex order, geometric parameters, and rotation angles + coeffs_h(tri, 0) = static_cast(order[0]); + coeffs_h(tri, 1) = static_cast(order[1]); + coeffs_h(tri, 2) = static_cast(order[2]); + coeffs_h(tri, 3) = a; + coeffs_h(tri, 4) = b; + coeffs_h(tri, 5) = c; + coeffs_h(tri, 6) = sin_theta; + coeffs_h(tri, 7) = cos_theta; + + // Compute coefficients + Real coeffs_tri[18][20]; + computeReducedQuinticCoefficients(a, b, c, coeffs_tri); + + // Store in flattened format after order, geometric parameters, and rotation angles + for (int i = 0; i < 18; i++) { + for (int j = 0; j < 20; j++) { + coeffs_h(tri, 8 + i * 20 + j) = coeffs_tri[i][j]; + } + } + } + + // Copy from host mirror to device + Kokkos::deep_copy(coeffs_d, coeffs_h); + + return coeffs_d; +} +} // namespace MeshField + +#endif // MESHFIELD_REDUCEDQUINTIC_HPP diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index 79b6b3f..c128d2c 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -225,5 +225,264 @@ struct QuadraticTetrahedronShape { } }; +/** + * @brief Helper functions for reduced quintic coordinate transformations + */ +namespace ReducedQuinticHelpers { + /** + * @brief Transform barycentric coordinates to reduced quintic local coordinates + * + * local coordinate system (origin at foot of perpendicular): + * - v0 is at (-b, 0) + * - v1 is at (a, 0) + * - v2 is at (0, c) + * + * @param xi Barycentric coordinates [λ0, λ1, λ2] + * @param a Distance from origin to v1 + * @param b Distance from origin to v0 + * @param c Perpendicular distance from origin to v2 + * @return Vector2 containing [xi_local, eta_local] + */ + KOKKOS_INLINE_FUNCTION + Vector2 barycentricToLocal( + Vector3 const& xi, + int const order[3], + Real a, + Real b, + Real c) + { + // Reorder barycentric coordinates into vertex ordering + const Real lambda0 = xi[order[0]]; + const Real lambda1 = xi[order[1]]; + const Real lambda2 = xi[order[2]]; + + const Real xi_local = a * lambda1 - b * lambda0; + const Real eta_local = c * lambda2; + + return {xi_local, eta_local}; + } + + /** + * @brief Get polynomial index for reduced quintic basis + * + * Returns [i,j] for xi^i * eta^j terms where i+j <= 5 + * Total of 21 possible terms, but we use 20 (the 21st is constrained) + * + * @param idx Index from 0 to 19 + * @return Pair of integers [xi_power, eta_power] + */ + KOKKOS_INLINE_FUNCTION + constexpr Kokkos::Array getReducedQuinticPolyIdx(int idx) { + // Order: (0,0), (1,0), (0,1), (2,0), (1,1), (0,2), ... + // Pattern: for each total degree d from 0 to 5, enumerate (i,j) where i+j=d + constexpr Kokkos::Array indices[20] = { + {0,0}, {1,0}, {0,1}, {2,0}, {1,1}, {0,2}, {3,0}, {2,1}, {1,2}, {0,3}, + {4,0}, {3,1}, {2,2}, {1,3}, {0,4}, {5,0}, {3,2}, {2,3}, {1,4}, {0,5} + }; + return indices[idx]; + } +} // namespace ReducedQuinticHelpers + +/** + * @brief Reduced quintic triangle element + * + * This element uses 18 nodes (3 vertices × 6 DOFs per vertex): + * - DOFs per vertex: [value, ∂/∂x, ∂/∂y, ∂²/∂x², ∂²/∂x∂y, ∂²/∂y²] + * - Polynomial order: 5 (20-term basis: xi^i * eta^j for i+j ≤ 5) + * + * COORDINATE TRANSFORMATION: + * This element uses element-specific local coordinates based on triangle geometry. + * + * The origin is placed at the foot of perpendicular from v2 onto the v0-v1 edge. + * - v0 is at (-b, 0) in local coords, where b = distance from origin to v0 + * - v1 is at (a, 0) in local coords, where a = distance from origin to v1 + * - v2 is at (0, c) in local coords, where c = perpendicular distance to v2 + * + * meshFields API uses barycentric coordinates: (λ0, λ1, λ2) where λ0+λ1+λ2=1 + * Passed as xi[0]=λ0, xi[1]=λ1, xi[2]=λ2. + * + * Transformation: + * xi_local = a*λ1 - b*λ0 + * eta_local = c*λ2 + * + * Applied automatically via helper functions: + * barycentricToLocal() and localToBarycentricGradient() + * + * The geometric parameters (a, b, c) are stored with the coefficients and + * retrieved during evaluation. Shape function coefficients are computed by + * solving a 20×20 linear system based on boundary conditions. + * + */ +struct ReducedQuinticTriangleShape { + static const size_t numNodes = 18; // 3 vertices × 6 DOFs per vertex + static const size_t meshEntDim = 2; + constexpr static Mesh_Topology DofHolders[1] = {Vertex}; + constexpr static size_t NumDofHolders[1] = {3}; // 3 vertices + constexpr static size_t DofsPerHolder[1] = {6}; // 6 DOFs per vertex + constexpr static size_t Order = 5; + + KOKKOS_INLINE_FUNCTION + Kokkos::Array getValues(Vector3 const &xi, + Kokkos::View elemCoeffs) const { + assert(greaterThanOrEqualZero(xi)); + assert(sumsToOne(xi)); + + // Extract geometric parameters from coefficient array + // elemCoeffs layout: [order[0], order[1], order[2], a, b, c, sin_theta, cos_theta, coeff_0_0, ..., coeff_17_19] + const int order[3] = {static_cast(elemCoeffs(0)), + static_cast(elemCoeffs(1)), + static_cast(elemCoeffs(2))}; + const Real a = elemCoeffs(3); + const Real b = elemCoeffs(4); + const Real c = elemCoeffs(5); + + // Transform barycentric to local coordinates + const auto local = ReducedQuinticHelpers::barycentricToLocal(xi, order, a, b, c); + const Real xi_local = local[0]; + const Real eta_local = local[1]; + + // Compute polynomial basis: xi_local^i * eta_local^j + Real xi_pow[6], eta_pow[6]; + xi_pow[0] = 1.0; eta_pow[0] = 1.0; + for (int i = 1; i < 6; i++) { + xi_pow[i] = xi_pow[i-1] * xi_local; + eta_pow[i] = eta_pow[i-1] * eta_local; + } + + // Evaluate shape functions using precomputed coefficients + Kokkos::Array N_reordered; + + for (size_t k = 0; k < numNodes; k++) { + N_reordered[k] = 0.0; + + for (int i = 0; i < 20; i++) { + const auto poly = ReducedQuinticHelpers::getReducedQuinticPolyIdx(i); + const int xi_idx = poly[0]; + const int eta_idx = poly[1]; + + N_reordered[k] += + elemCoeffs(8 + k * 20 + i) * + xi_pow[xi_idx] * + eta_pow[eta_idx]; + } + } + + // Convert back to meshFields vertex ordering + Kokkos::Array N; + + for (int v = 0; v < 3; ++v) { + const int orig_v = order[v]; + + for (int d = 0; d < 6; ++d) { + N[orig_v * 6 + d] = + N_reordered[v * 6 + d]; + } + } + + return N; + } + + KOKKOS_INLINE_FUNCTION + Kokkos::Array getLocalGradients(Vector3 const &xi, + Kokkos::View elemCoeffs) const { + assert(greaterThanOrEqualZero(xi)); + assert(sumsToOne(xi)); + + // Extract geometric parameters + const int order[3] = {static_cast(elemCoeffs(0)), + static_cast(elemCoeffs(1)), + static_cast(elemCoeffs(2))}; + const Real a = elemCoeffs(3); + const Real b = elemCoeffs(4); + const Real c = elemCoeffs(5); + + // Transform barycentric to local coordinates + const auto local = ReducedQuinticHelpers::barycentricToLocal(xi, order, a, b, c); + const Real xi_local = local[0]; + const Real eta_local = local[1]; + + // Compute polynomial basis and derivatives in local coordinates + Real xi_pow[6], eta_pow[6]; + Real dxi_pow[6], deta_pow[6]; + + xi_pow[0] = 1.0; eta_pow[0] = 1.0; + dxi_pow[0] = 0.0; deta_pow[0] = 0.0; + + for (int i = 1; i < 6; i++) { + xi_pow[i] = xi_pow[i-1] * xi_local; + eta_pow[i] = eta_pow[i-1] * eta_local; + dxi_pow[i] = i * xi_pow[i-1]; + deta_pow[i] = i * eta_pow[i-1]; + } + + // Compute matrix for local to barycentric gradient chain rule + Real J[2][2] = {{0.0, 0.0}, {0.0, 0.0}}; + for (int col = 0; col < 2; col++) { + // d(lambda_k)/d(xi[col]) for k = 0, 1, 2 + Real dlambda[3]; + for (int k = 0; k < 3; k++) { + if (order[k] == col) dlambda[k] = 1.0; + else if (order[k] == 2) dlambda[k] = -1.0; + else dlambda[k] = 0.0; + } + + // d(xi_local)/d(xi[col]) = a*dlambda[1] - b*dlambda[0] + // d(eta_local)/d(xi[col]) = c*dlambda[2] + J[0][col] = a * dlambda[1] - b * dlambda[0]; + J[1][col] = c * dlambda[2]; + } + + // Evaluate gradients + Kokkos::Array grad_reordered; + + for (size_t k = 0; k < numNodes; k++) { + Real dN_dxi_local = 0.0; + Real dN_deta_local = 0.0; + + for (int i = 0; i < 20; i++) { + const auto poly = ReducedQuinticHelpers::getReducedQuinticPolyIdx(i); + const int xi_idx = poly[0]; + const int eta_idx = poly[1]; + const Real coeff = elemCoeffs(8 + k * 20 + i); + + if (xi_idx > 0) + dN_dxi_local += + coeff * + dxi_pow[xi_idx] * + eta_pow[eta_idx]; + + if (eta_idx > 0) + dN_deta_local += + coeff * + xi_pow[xi_idx] * + deta_pow[eta_idx]; + } + + // Apply chain rule to transform local gradients to barycentric gradients + grad_reordered[k][0] = + dN_dxi_local * J[0][0] + + dN_deta_local * J[1][0]; + + grad_reordered[k][1] = + dN_dxi_local * J[0][1] + + dN_deta_local * J[1][1]; + } + + // Convert back to meshFields vertex ordering + Kokkos::Array gradN_bary; + + for (int v = 0; v < 3; ++v) { + const int orig_v = order[v]; + + for (int d = 0; d < 6; ++d) { + gradN_bary[orig_v * 6 + d] = + grad_reordered[v * 6 + d]; + } + } + + return gradN_bary; + } +}; + } // namespace MeshField #endif diff --git a/test/testOmegahTriReducedQuintic.cpp b/test/testOmegahTriReducedQuintic.cpp new file mode 100644 index 0000000..99e0aec --- /dev/null +++ b/test/testOmegahTriReducedQuintic.cpp @@ -0,0 +1,305 @@ +#include "KokkosController.hpp" +#include "MeshField.hpp" +#include "MeshField_Element.hpp" +#include "MeshField_Fail.hpp" +#include "MeshField_For.hpp" +#include "MeshField_ShapeField.hpp" +#include "Omega_h_build.hpp" +#include "Omega_h_file.hpp" +#include "Omega_h_simplex.hpp" +#include +#include +#include + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; + +struct LinearFunction { + KOKKOS_INLINE_FUNCTION + MeshField::Real operator()(MeshField::Real x, MeshField::Real y) const { + return 2.0 * x + y; + } + // ∂f/∂x = 2.0, ∂f/∂y = 1.0 + // ∂²f/∂x² = 0.0, ∂²f/∂xy = 0.0, ∂²f/∂y² = 0.0 + static constexpr MeshField::Real dfdx = 2.0; + static constexpr MeshField::Real dfdy = 1.0; + static constexpr MeshField::Real d2fdx2 = 0.0; + static constexpr MeshField::Real d2fdxy = 0.0; + static constexpr MeshField::Real d2fdy2 = 0.0; +}; + +struct QuadraticFunction { + KOKKOS_INLINE_FUNCTION + MeshField::Real operator()(MeshField::Real x, MeshField::Real y) const { + return (x * x) + (2.0 * y); + } + // ∂f/∂x = 2x, ∂f/∂y = 2.0 + // ∂²f/∂x² = 2.0, ∂²f/∂xy = 0.0, ∂²f/∂y² = 0.0 +}; + +Omega_h::Mesh createMeshTri18(Omega_h::Library &lib) { + auto world = lib.world(); + const auto family = OMEGA_H_SIMPLEX; + auto len = 1.0; + return Omega_h::build_box(world, family, len, len, 0.0, 3, 3, 0); +} + +struct TestCoords { + Kokkos::View coords; + size_t NumPtsPerElem; + std::string name; +}; + +template +bool checkResult(Omega_h::Mesh &mesh, Result &result, CoordField coordField, + TestCoords testCase, AnalyticFunction func, size_t numComp) { + // Create explicit local copy to ensure proper capture in KOKKOS_LAMBDA + const size_t npts_per_elem = testCase.NumPtsPerElem; + const size_t num_comp = numComp; + + MeshField::FieldElement fcoords( + mesh.nfaces(), coordField, MeshField::LinearTriangleCoordinateShape(), + MeshField::Omegah::LinearTriangleToVertexField(mesh)); + auto globalCoords = + MeshField::evaluate(fcoords, testCase.coords, npts_per_elem); + + // Host-side debugging output + std::cout << "\n=== DEBUG checkResult for " << testCase.name << " ===" << std::endl; + std::cout << " numPtsPerElem = " << npts_per_elem << std::endl; + std::cout << " numComponents = " << num_comp << std::endl; + std::cout << " mesh.nfaces() = " << mesh.nfaces() << std::endl; + std::cout << " result.extent(0) = " << result.extent(0) << std::endl; + std::cout << " result.extent(1) = " << result.extent(1) << std::endl; + std::cout << " globalCoords.extent(0) = " << globalCoords.extent(0) << std::endl; + std::cout << " Expected total points = " << mesh.nfaces() * npts_per_elem << std::endl; + std::cout << " MachinePrecision = " << MeshField::MachinePrecision << std::endl; + + MeshField::LO numErrors = 0; + Kokkos::parallel_reduce( + "checkResult", mesh.nfaces(), + KOKKOS_LAMBDA(const int &ent, MeshField::LO &lerrors) { + const auto first = ent * npts_per_elem; + const auto last = first + npts_per_elem; + + // Bounds checking + if (first >= result.extent(0) || last > result.extent(0)) { + Kokkos::printf("ERROR: ent=%d, first=%d, last=%d exceeds result.extent(0)=%d\n", + ent, (int)first, (int)last, (int)result.extent(0)); + lerrors += 1; + return; + } + + if (first >= globalCoords.extent(0) || last > globalCoords.extent(0)) { + Kokkos::printf("ERROR: ent=%d, first=%d, last=%d exceeds globalCoords.extent(0)=%d\n", + ent, (int)first, (int)last, (int)globalCoords.extent(0)); + lerrors += 1; + return; + } + + for (auto pt = first; pt < last; pt++) { + const auto x = globalCoords(pt, 0); + const auto y = globalCoords(pt, 1); + const auto expected = func(x, y); + + // Debug output for first few elements + if (ent < 3 && pt == first) { + Kokkos::printf("Element %d: pt=%d, x=%.6f, y=%.6f, expected=%.6f, result(pt,0)=%.6f\n", + ent, (int)pt, x, y, expected, result(pt, 0)); + } + + for (size_t i = 0; i < num_comp; ++i) { + const auto computed = result(pt, i); + const auto diff = Kokkos::fabs(computed - expected); + MeshField::LO isError = 0; + if (diff > MeshField::MachinePrecision) { + isError = 1; + Kokkos::printf( + "result for elm %d, pt %d (local %d), comp %d: expected " + "%.15f (x=%.6f, y=%.6f) computed %.15f, diff=%.15e, tol=%.15e\n", + ent, (int)pt, (int)(pt - first), (int)i, expected, x, y, computed, diff, MeshField::MachinePrecision); + } + lerrors += isError; + } + } + }, + numErrors); + + std::cout << " Total errors found: " << numErrors << std::endl; + return (numErrors > 0); +} + +template +void setReducedQuinticDOFs(Omega_h::Mesh &mesh, AnalyticFunction func, + ShapeField field) { + const auto MeshDim = mesh.dim(); + auto coords = mesh.coords(); + auto setFieldAtVertices = KOKKOS_LAMBDA(const int &vtx) { + const auto x = coords[vtx * MeshDim]; + const auto y = coords[vtx * MeshDim + 1]; + // value: f(x,y) + field(vtx, 0, 0, MeshField::Vertex) = func(x, y); + // first derivatives + if constexpr (std::is_same_v) { + field(vtx, 0, 1, MeshField::Vertex) = AnalyticFunction::dfdx; + field(vtx, 0, 2, MeshField::Vertex) = AnalyticFunction::dfdy; + field(vtx, 0, 3, MeshField::Vertex) = AnalyticFunction::d2fdx2; + field(vtx, 0, 4, MeshField::Vertex) = AnalyticFunction::d2fdxy; + field(vtx, 0, 5, MeshField::Vertex) = AnalyticFunction::d2fdy2; + } else { + // QuadraticFunction: f(x,y) = x^2 + 2y + // ∂f/∂x = 2x, ∂f/∂y = 2 + // ∂²f/∂x² = 2, ∂²f/∂xy = 0, ∂²f/∂y² = 0 + field(vtx, 0, 1, MeshField::Vertex) = 2.0 * x; + field(vtx, 0, 2, MeshField::Vertex) = 2.0; + field(vtx, 0, 3, MeshField::Vertex) = 2.0; + field(vtx, 0, 4, MeshField::Vertex) = 0.0; + field(vtx, 0, 5, MeshField::Vertex) = 0.0; + } + }; + MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()}, + setFieldAtVertices, "setReducedQuinticDOFs"); +} + +template +Kokkos::View +createElmAreaCoords(size_t numElements, + Kokkos::Array coords) { + Kokkos::View lc("localCoords", + numElements * NumPtsPerElem); + Kokkos::parallel_for( + "setLocalCoords", numElements, KOKKOS_LAMBDA(const int &elm) { + for (size_t pt = 0; pt < NumPtsPerElem; pt++) { + lc(elm * NumPtsPerElem + pt, 0) = coords[pt * 3 + 0]; + lc(elm * NumPtsPerElem + pt, 1) = coords[pt * 3 + 1]; + lc(elm * NumPtsPerElem + pt, 2) = coords[pt * 3 + 2]; + } + }); + return lc; +} + +void doFail(std::string_view order, std::string_view function, + std::string_view location, std::string_view numComp) { + std::stringstream ss; + ss << order << " field evaluation with " << numComp << " components and " + << function << " analytic function at " << location << " points failed\n"; + std::string msg = ss.str(); + MeshField::fail(msg); +} + +template typename Controller, + typename TestCaseType, typename FunctionType> +void runTest(Omega_h::Mesh &mesh, + MeshField::OmegahMeshField &omf, + TestCaseType testCase, FunctionType function) { + using functionType = decltype(function); + using ViewType = decltype(testCase.coords); + auto field = omf.template CreateLagrangeField(); + using FieldType = decltype(field); + setReducedQuinticDOFs(mesh, function, field); + // verify field + Kokkos::parallel_for( + "printField", mesh.nverts(), + KOKKOS_LAMBDA(const int &vtx) { + Kokkos::printf("vtx %d: value %f, dfdx %f, dfdy %f, d2fdx2 %f, d2fdxy %f, " + "d2fdy2 %f\n", + vtx, field(vtx, 0, 0, MeshField::Vertex), + field(vtx, 0, 1, MeshField::Vertex), + field(vtx, 0, 2, MeshField::Vertex), + field(vtx, 0, 3, MeshField::Vertex), + field(vtx, 0, 4, MeshField::Vertex), + field(vtx, 0, 5, MeshField::Vertex)); + }); + auto result = omf.template triangleReducedQuinticEval( + testCase.coords, testCase.NumPtsPerElem, field); + auto failed = checkResult(mesh, result, omf.getCoordField(), testCase, + decltype(function){}, numComponents); + if (failed) { + std::string functionErr; + if constexpr (std::is_same_v) { + functionErr = "linear"; + } else { + functionErr = "quadratic"; + } + doFail("reducedQuintic", functionErr, testCase.name, std::to_string(numComponents)); + } +} + +template