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..d50fe48 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 { + 0, // node within vertex dof holder + vertexDof, // field component (the DOF type) + vtx, + MeshField::Vertex + }; + } +}; + template auto getTriangleElement(Omega_h::Mesh &mesh) { static_assert(ShapeOrder == 1 || ShapeOrder == 2); if constexpr (ShapeOrder == 1) { @@ -256,6 +317,68 @@ template auto getTriangleElement(Omega_h::Mesh &mesh) { QuadraticTriangleToField(mesh)}; } } + +// ReducedQuintic triangle element with precomputed coefficients +// Functor for extracting triangle vertex coordinates on device +struct ExtractTriCoords { + Omega_h::LOs triVerts_d; + Omega_h::Reals coords_d; + Kokkos::View triCoords_d; + + ExtractTriCoords(Omega_h::LOs triVerts_d_, Omega_h::Reals coords_d_, + Kokkos::View triCoords_d_) + : triVerts_d(triVerts_d_), coords_d(coords_d_), triCoords_d(triCoords_d_) {} + + KOKKOS_INLINE_FUNCTION + void operator()(int tri) const { + 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_d[tri * triToVtxDegree + localVtxIdx]; + + triCoords_d(tri, vi * 2 + 0) = coords_d[vtx * 2 + 0]; // x + triCoords_d(tri, vi * 2 + 1) = coords_d[vtx * 2 + 1]; // y + } + } +}; + +// Result type for getReducedQuinticTriangleElement +struct ReducedQuinticTriangleElementResult { + MeshField::ReducedQuinticTriangleShape shp; + Omegah::ReducedQuinticTriangleToField map; + Kokkos::View coeffs; +}; + +inline ReducedQuinticTriangleElementResult +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(); + + // Allocate device view for triangle vertex coordinates [numTri][6] + // Layout: tri*6 + v*2 + d where v in {0,1,2} and d in {0,1} for (x,y) + Kokkos::View triCoords_d("triCoords_device", numTri, 6); + + // Extract triangle vertex coordinates entirely on device using functor + Kokkos::parallel_for( + "extractTriCoords", numTri, + ExtractTriCoords(triVerts_d, coords_d, triCoords_d)); + + // Precompute coefficients for all triangles entirely on device + auto elemCoeffs = MeshField::precomputeReducedQuinticCoefficients(numTri, triCoords_d); + + return ReducedQuinticTriangleElementResult{ + MeshField::ReducedQuinticTriangleShape(), + Omegah::ReducedQuinticTriangleToField(mesh), + elemCoeffs}; +} template auto getTetrahedronElement(Omega_h::Mesh &mesh) { static_assert(ShapeOrder == 1 || ShapeOrder == 2); if constexpr (ShapeOrder == 1) { @@ -360,6 +483,38 @@ class OmegahMeshField { return eval; } + // Evaluate a ReducedQuintic field at the specified local coordinates for each triangle. + // + // NOTE: This is separate from triangleLocalPointEval because the ReducedQuintic + // element has a different DOF layout and shape function. triangleLocalPointEval + // handles standard elements, while ReducedQuintic uses + // a different set of degrees of freedom and shape functions. + 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..72a1e7f 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]; + } } } } @@ -207,6 +273,28 @@ struct FieldElement { return c; } + // Geometry node array: uses GeometryShape::numNodes + // For isoparametric elements, this equals ShapeType::numNodes + // For non-isoparametric (e.g., ReducedQuintic), this equals the geometry node count + // This would replace the getNodeValues function for Jacobian computation + using GeomNodeArray = + Kokkos::Array::type, + ShapeType::meshEntDim * ShapeType::GeometryShape::numNodes>; + KOKKOS_INLINE_FUNCTION GeomNodeArray getGeometryNodeValues(int ent) const { + GeomNodeArray c; + for (auto topo : elm2dof.getTopology()) { // element topology + for (size_t gni = 0; gni < ShapeType::GeometryShape::numNodes; ++gni) { + for (size_t d = 0; d < ShapeType::meshEntDim; ++d) { + auto map = elm2dof(gni, d, ent, topo); + const auto fval = + field(map.entity, map.node, map.component, map.topo); + c[gni * ShapeType::meshEntDim + d] = fval; + } + } + } + return c; + } + /** * @brief * compute the Jacobian of an edge @@ -222,10 +310,11 @@ struct FieldElement { KOKKOS_INLINE_FUNCTION Real getJacobian1d(int ent) const { assert(ent >= 0); assert(static_cast(ent) < numMeshEnts); - const auto nodalGradients = shapeFn.getLocalGradients(); - const auto nodeValues = getNodeValues(ent); + typename ShapeType::GeometryShape geomShape{}; + const auto nodalGradients = geomShape.getLocalGradients(); + const auto nodeValues = getGeometryNodeValues(ent); auto g = nodalGradients[0] * nodeValues[0]; - for (size_t i = 1; i < shapeFn.numNodes; ++i) { + for (size_t i = 1; i < ShapeType::GeometryShape::numNodes; ++i) { g = g + nodalGradients[i] * nodeValues[i]; } return g; @@ -383,21 +472,24 @@ struct FieldElement { Kokkos::View res("result", numPts, MeshEntDim, MeshEntDim); Kokkos::deep_copy(res, 0.0); // initialize all entries to zero - // fill the views of node coordinates and node gradients - Kokkos::View nodeCoords( + // Use GeometryShape for the geometry mapping gradient (supports non-isoparametric elements) + typename ShapeType::GeometryShape geomShape{}; + const auto geomGrad = geomShape.getLocalGradients(); + + // fill the views of geometry node coordinates and geometry node gradients + Kokkos::View nodeCoords( "nodeCoords", numPts); - Kokkos::View nodalGradients( + Kokkos::View nodalGradients( "nodalGradients", numPts); - const auto grad = shapeFn.getLocalGradients(); Kokkos::parallel_for( numMeshEnts, KOKKOS_CLASS_LAMBDA(const int ent) { - const auto vals = getNodeValues(ent); - assert(vals.size() == MeshEntDim * ShapeType::numNodes); + const auto vals = getGeometryNodeValues(ent); + assert(vals.size() == MeshEntDim * ShapeType::GeometryShape::numNodes); for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) { - for (size_t node = 0; node < ShapeType::numNodes; node++) { + for (size_t node = 0; node < ShapeType::GeometryShape::numNodes; node++) { for (size_t d = 0; d < MeshEntDim; d++) { nodeCoords(pt, node, d) = vals[node * MeshEntDim + d]; - nodalGradients(pt, node, d) = grad[node * MeshEntDim + d]; + nodalGradients(pt, node, d) = geomGrad[node * MeshEntDim + d]; } } } @@ -408,7 +500,7 @@ struct FieldElement { // TODO use nested parallel for? for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) { auto A = Kokkos::subview(res, pt, Kokkos::ALL(), Kokkos::ALL()); - for (size_t node = 0; node < ShapeType::numNodes; node++) { + for (size_t node = 0; node < ShapeType::GeometryShape::numNodes; node++) { auto a = Kokkos::subview(nodalGradients, pt, node, Kokkos::ALL()); auto b = Kokkos::subview(nodeCoords, pt, node, Kokkos::ALL()); diff --git a/src/MeshField_ReducedQuintic.hpp b/src/MeshField_ReducedQuintic.hpp new file mode 100644 index 0000000..8251b24 --- /dev/null +++ b/src/MeshField_ReducedQuintic.hpp @@ -0,0 +1,430 @@ +#ifndef MESHFIELD_REDUCEDQUINTIC_HPP +#define MESHFIELD_REDUCEDQUINTIC_HPP + +#include +#include +#include +#include + +namespace MeshField { + +/** + * @brief Device-compatible 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. + * + * @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 + */ +template +KOKKOS_INLINE_FUNCTION +int solveLU_internal(int n, int nrhs, Real* A, int lda, Real* B, int ldb) { + // Use fixed-size array for pivoting (n=20 maximum, which fits all our use cases) + // We use a stack array since this is device-compatible + int pivot[20]; + const Real eps = 1e-14; + + // LU decomposition with partial pivoting + for (int k = 0; k < n; k++) { + // Find pivot + int maxRow = k; + Real maxVal = Kokkos::fabs(A[k * lda + k]); + for (int i = k + 1; i < n; i++) { + Real val = Kokkos::fabs(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++) { + Real tmp = A[k * lda + j]; + A[k * lda + j] = A[maxRow * lda + j]; + A[maxRow * lda + j] = tmp; + } + } + + // 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++) { + Real tmp = B[k * ldb + rhs]; + B[k * ldb + rhs] = B[pivot[k] * ldb + rhs]; + B[pivot[k] * ldb + rhs] = tmp; + } + } + } + + // 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 Device-compatible function to 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 + */ +template +KOKKOS_INLINE_FUNCTION +void rotateDof(Real dofs_p[6], Real sin_theta_p, Real cos_theta_p) +{ + const Real s = sin_theta_p; + const Real c = cos_theta_p; + const Real ss = s * s; + const Real cc = c * c; + const Real sc = s * c; + + const Real 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_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 Device-compatible function to reorder triangle vertices to put longest edge along local xi-axis + * + * @param coords Triangle vertex coordinates [3][2] + * @param order Output vertex order [3] - order[i] gives original index of reordered vertex i + */ +template +KOKKOS_INLINE_FUNCTION +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 Device-compatible function to compute geometric parameters for reduced quintic triangle element + * + * @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 + */ +template +KOKKOS_INLINE_FUNCTION +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]; + } + + // 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 = Kokkos::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 = Kokkos::sqrt(vec_a[0]*vec_a[0] + vec_a[1]*vec_a[1]); // distance to v0 + a = Kokkos::sqrt(vec_b[0]*vec_b[0] + vec_b[1]*vec_b[1]); // distance to v1 + c = Kokkos::sqrt(vec_c[0]*vec_c[0] + vec_c[1]*vec_c[1]); // distance to v2 + + // Angle of v0-v1 edge + sin_theta = ab[1] / ablen; + cos_theta = ab[0] / ablen; +} + +/** + * @brief Device-compatible function to compute coefficient matrix for reduced quintic element + * + * @param a, b, c Geometric parameters from computeReducedQuinticGeometry + * @param coeffs Output: coefficient matrix [18][20] + */ +template +KOKKOS_INLINE_FUNCTION +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) + 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); + + // 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, entirely on device + * + * @param numTriangles Number of triangles in the mesh + * @param triCoords_d Triangle vertex coordinates on device, shape [numTriangles][3][2] (flattened: tri*6 + v*2 + d) + * @return Kokkos::View with coefficients, shape [numTriangles][8 + 18*20] + * Each row: [order[0], order[1], order[2], a, b, c, sin_theta, cos_theta, coeff_0_0, ..., coeff_17_19] + */ +inline Kokkos::View precomputeReducedQuinticCoefficients( + int numTriangles, + Kokkos::View triCoords_d) +{ + // 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); + + // Compute for each triangle entirely on device + Kokkos::parallel_for( + "precomputeReducedQuinticCoefficients", numTriangles, + KOKKOS_LAMBDA(int tri) { + Real coords[3][2]; + for (int v = 0; v < 3; v++) { + coords[v][0] = triCoords_d(tri, v * 2 + 0); + coords[v][1] = triCoords_d(tri, 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_d(tri, 0) = static_cast(order[0]); + coeffs_d(tri, 1) = static_cast(order[1]); + coeffs_d(tri, 2) = static_cast(order[2]); + coeffs_d(tri, 3) = a; + coeffs_d(tri, 4) = b; + coeffs_d(tri, 5) = c; + coeffs_d(tri, 6) = sin_theta; + coeffs_d(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_d(tri, 8 + i * 20 + j) = coeffs_tri[i][j]; + } + } + }); + + return coeffs_d; +} +} // namespace MeshField + +#endif // MESHFIELD_REDUCEDQUINTIC_HPP \ No newline at end of file diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index 79b6b3f..1e82e19 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -43,6 +43,7 @@ struct LinearEdgeShape { static const size_t meshEntDim = 1; constexpr static Mesh_Topology DofHolders[1] = {Vertex}; constexpr static size_t Order = 1; + using GeometryShape = LinearEdgeShape; KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector2 const &xi) const { @@ -69,6 +70,7 @@ struct LinearTriangleShape { static const size_t meshEntDim = 2; constexpr static Mesh_Topology DofHolders[1] = {Vertex}; constexpr static size_t Order = 1; + using GeometryShape = LinearTriangleShape; KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector3 const &xi) const { @@ -96,6 +98,7 @@ struct LinearTriangleCoordinateShape { static const size_t meshEntDim = 2; constexpr static Mesh_Topology DofHolders[1] = {Vertex}; constexpr static size_t Order = 1; + using GeometryShape = LinearTriangleCoordinateShape; KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector3 const &xi) const { @@ -114,6 +117,7 @@ struct LinearTetrahedronShape { static const size_t meshEntDim = 3; constexpr static Mesh_Topology DofHolders[1] = {Vertex}; constexpr static size_t Order = 1; + using GeometryShape = LinearTetrahedronShape; KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector4 const &xi) const { @@ -144,6 +148,7 @@ struct QuadraticTriangleShape { constexpr static size_t NumDofHolders[2] = {3, 3}; constexpr static size_t DofsPerHolder[2] = {1, 1}; constexpr static size_t Order = 2; + using GeometryShape = QuadraticTriangleShape; KOKKOS_INLINE_FUNCTION Kokkos::Array getValues(Vector3 const &xi) const { @@ -183,6 +188,7 @@ struct QuadraticTetrahedronShape { constexpr static size_t NumDofHolders[2] = {4, 6}; constexpr static size_t DofsPerHolder[2] = {1, 1}; constexpr static size_t Order = 2; + using GeometryShape = QuadraticTetrahedronShape; // ordering taken from mfem // see mfem/mfem fem/fe/fe_fixed_order.cpp @597cba8 KOKKOS_INLINE_FUNCTION @@ -225,5 +231,283 @@ 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; + using GeometryShape = LinearTriangleShape; + + 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); + + // Reorder barycentric from meshFields pipeline convention to coords-index convention + const Real xi0_ordered = 1.0 - xi[0] - xi[1]; // coords[0] + const Real xi1_ordered = xi[0]; // coords[1] + const Real xi2_ordered = xi[1]; // coords[2] + Vector3 xi_ordered = {xi0_ordered, xi1_ordered, xi2_ordered}; + + // Transform barycentric to local coordinates + const auto local = ReducedQuinticHelpers::barycentricToLocal(xi_ordered, 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; + } + + // Computes the local gradients of the Reduced Quintic shape functions. + // + // These gradients are not currently used in the MeshField evaluation pipeline + // because the Reduced Quintic element uses a linear geometric mapping. The + // Jacobian is therefore computed from the linear geometry rather than from the + // Reduced Quintic shape functions. + 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); + + // Reorder barycentric from meshFields pipeline convention to coords-index convention + const Real xi0_ordered = 1.0 - xi[0] - xi[1]; // coords[0] + const Real xi1_ordered = xi[0]; // coords[1] + const Real xi2_ordered = xi[1]; // coords[2] + Vector3 xi_ordered = {xi0_ordered, xi1_ordered, xi2_ordered}; + + // Transform barycentric to local coordinates + const auto local = ReducedQuinticHelpers::barycentricToLocal(xi_ordered, 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..db73a13 --- /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 > 10 * 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