From a3e920801637d1ab807693d1d03a1e233ebac46f Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Tue, 6 Oct 2026 02:43:06 -0400 Subject: [PATCH 01/19] add function to perform interpolation onto dofs This performs interpolation to target elements. We use ownership to avoid having multiple element write to the same dofs on shared entities. This interpolation capability will be used to support generic coordinate field initialization for higher order elements. --- CMakeLists.txt | 2 + src/MeshField_Element.hpp | 84 ++++++++++++++++ src/MeshField_ShapeField.hpp | 41 ++++++-- test/testInterpolate.cpp | 179 +++++++++++++++++++++++++++++++++++ 4 files changed, 297 insertions(+), 9 deletions(-) create mode 100644 test/testInterpolate.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 45dc805..f295708 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(Interpolate test/testInterpolate.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(Interpolate ./Interpolate) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index cb2890d..b729b24 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -599,5 +599,89 @@ Kokkos::View getJacobianDeterminants(FieldElement &fes, Matrices J) { return fes.getJacobianDeterminants(J); } +/** + * @brief + * For each dof holder of a field, the element + * that owns it; see getDofHolderOwningEntities + */ +struct DofHolderOwningEntities { + static constexpr size_t NumTopologies = Pyramid + 1; + Kokkos::Array offset; + Kokkos::View owner; + + /** @brief the index of the element that owns the given dof holder */ + KOKKOS_INLINE_FUNCTION LO operator()(Mesh_Topology topo, LO entity) const { + return owner(offset[topo] + entity); + } +}; + +/** + * @brief + * the owning entity of a dof that is shared, i.e., vertex or edge nodes + * is the shared entity with the lowest id + */ +template +DofHolderOwningEntities getDofHolderOwningEntities(FieldElement const &fes) { + constexpr size_t numNodes = decltype(fes.shapeFn)::numNodes; + constexpr size_t numTopologies = DofHolderOwningEntities::NumTopologies; + const auto &meshInfo = fes.field.meshInfo; + const LO numElements = fes.numMeshEnts; + DofHolderOwningEntities owners; + owners.offset[0] = 0; + for (size_t t = 0; t < numTopologies; ++t) + owners.offset[t + 1] = + owners.offset[t] + meshInfo.numEntities(static_cast(t)); + owners.owner = Kokkos::View( + Kokkos::view_alloc(Kokkos::WithoutInitializing, "dofHolderOwner"), + owners.offset[numTopologies]); + Kokkos::deep_copy(owners.owner, numElements); + const auto offset = owners.offset; + const auto owner = owners.owner; + Kokkos::parallel_for( + "getDofHolderOwningEntities", numElements, KOKKOS_LAMBDA(const int ent) { + for (auto topo : fes.elm2dof.getTopology()) { + for (size_t n = 0; n < numNodes; ++n) { + const auto map = fes.elm2dof(n, 0, ent, topo); + Kokkos::atomic_min(&owner(offset[map.topo] + map.entity), ent); + } + } + }); + return owners; +} + +/** + * @brief + * evaluates a source function at of the field element dof holder location + */ +template +void interpolate(FieldElement &target, Source const &source) { + constexpr size_t numNodes = decltype(target.shapeFn)::numNodes; + constexpr size_t meshEntDim = FieldElement::MeshEntDim; + constexpr size_t numComp = FieldElement::NumComponents; + // we get a list of owners so that we only perform evaluation once + // per shared dof + const auto owners = getDofHolderOwningEntities(target); + Kokkos::parallel_for( + "interpolate", target.numMeshEnts, KOKKOS_LAMBDA(const int ent) { + const auto nodeCoords = target.shapeFn.getNodeParametricCoords(); + for (auto topo : target.elm2dof.getTopology()) { + for (size_t n = 0; n < numNodes; ++n) { + const auto holder = target.elm2dof(n, 0, ent, topo); + if (owners(holder.topo, holder.entity) != ent) + continue; + Kokkos::Array xi; + for (size_t d = 0; d < meshEntDim; ++d) + xi[d] = nodeCoords[n * meshEntDim + d]; + const auto val = source(ent, xi); + for (size_t c = 0; c < numComp; ++c) { + const auto map = target.elm2dof(n, c, ent, topo); + target.field(map.entity, map.node, map.component, map.topo) = + val[c]; + } + } + } + }); +} + } // namespace MeshField #endif diff --git a/src/MeshField_ShapeField.hpp b/src/MeshField_ShapeField.hpp index 5175b02..6178a80 100644 --- a/src/MeshField_ShapeField.hpp +++ b/src/MeshField_ShapeField.hpp @@ -17,15 +17,38 @@ namespace MeshField { * On-process mesh metadata */ struct MeshInfo { - int numVtx; // entDim = 0 - int numEdge; // entDim = 1 - int numTri; // entDim = 2 - int numQuad; // entDim = 2 - int numTet; // entDim = 3 - int numHex; // entDim = 3 - int numPrism; // entDim = 3 - int numPyramid; // entDim = 3 - int dim; + int numVtx = 0; // entDim = 0 + int numEdge = 0; // entDim = 1 + int numTri = 0; // entDim = 2 + int numQuad = 0; // entDim = 2 + int numTet = 0; // entDim = 3 + int numHex = 0; // entDim = 3 + int numPrism = 0; // entDim = 3 + int numPyramid = 0; // entDim = 3 + int dim = 0; + + /** @brief runtime map of number of entities based on topology */ + int numEntities(Mesh_Topology topo) const { + switch (topo) { + case Vertex: + return numVtx; + case Edge: + return numEdge; + case Triangle: + return numTri; + case Quad: + return numQuad; + case Tetrahedron: + return numTet; + case Hexahedron: + return numHex; + case Prism: + return numPrism; + case Pyramid: + return numPyramid; + } + return 0; + } }; /** diff --git a/test/testInterpolate.cpp b/test/testInterpolate.cpp new file mode 100644 index 0000000..ffdcb94 --- /dev/null +++ b/test/testInterpolate.cpp @@ -0,0 +1,179 @@ +#include "KokkosController.hpp" +#include "MeshField.hpp" +#ifdef MESHFIELDS_ENABLE_CABANA +#include "CabanaController.hpp" +#endif +#include "Omega_h_build.hpp" +#include +#include +#include + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +using MeshField::Real; + +template auto getElement(Omega_h::Mesh &mesh) { + if constexpr (dim == 3) { + return MeshField::Omegah::getTetrahedronElement(mesh); + } else { + return MeshField::Omegah::getTriangleElement(mesh); + } +} + +template +KOKKOS_INLINE_FUNCTION Real quadratic(Kokkos::Array const &x) { + Real f = 1; + for (size_t d = 0; d < dim; ++d) + f += (d + 1) * x[d] + x[d] * x[d]; + return f; +} + +template +KOKKOS_INLINE_FUNCTION Real linear(Kokkos::Array const &x) { + Real f = 1; + for (size_t d = 0; d < dim; ++d) + f += (d + 1) * x[d]; + return f; +} + +/* + * Interpolate an analytic function, with numComp components + */ +template typename Controller> +bool testAnalytic( + Omega_h::Mesh &mesh, + MeshField::OmegahMeshField &omf) { + const auto [coordShp, coordMap] = getElement(mesh); + MeshField::FieldElement coordElm(mesh.nelems(), omf.getCoordField().field, + coordShp, coordMap); + auto fieldWithCtrlr = + omf.template CreateLagrangeField(); + const auto [shp, map] = getElement(mesh); + MeshField::FieldElement elm(mesh.nelems(), fieldWithCtrlr.field, shp, map); + + const auto f = KOKKOS_LAMBDA(Kokkos::Array const &x) { + return order == 1 ? linear(x) : quadratic(x); + }; + MeshField::interpolate( + elm, KOKKOS_LAMBDA(const int ent, Kokkos::Array const &xi) { + const auto x = coordElm.getValue(ent, xi); + Kokkos::Array val; + for (size_t c = 0; c < numComp; ++c) + val[c] = (c + 1) * f(x); + return val; + }); + + // points away from the nodes: the centroid and an off-center point + constexpr size_t numPts = 2; + Kokkos::View lc("localCoords", mesh.nelems() * numPts, dim); + Kokkos::parallel_for( + "setLocalCoords", mesh.nelems(), KOKKOS_LAMBDA(const int ent) { + for (size_t d = 0; d < dim; ++d) { + lc(ent * numPts, d) = 1.0 / (dim + 1); + lc(ent * numPts + 1, d) = 0.1 * (d + 1); + } + }); + const auto values = MeshField::evaluate(elm, lc, numPts); + const auto x = MeshField::evaluate(coordElm, lc, numPts); + MeshField::LO numErrors = 0; + Kokkos::parallel_reduce( + "checkAnalytic", values.extent(0), + KOKKOS_LAMBDA(const int pt, MeshField::LO &lerrors) { + Kokkos::Array xp; + for (size_t d = 0; d < dim; ++d) + xp[d] = x(pt, d); + for (size_t c = 0; c < numComp; ++c) + if (Kokkos::fabs(values(pt, c) - (c + 1) * f(xp)) > + MeshField::Epsilon) + ++lerrors; + }, + numErrors); + const bool pass = numErrors == 0; + std::cout << "testAnalytic(dim " << dim << ", order " << order << ", " + << numComp << " components): " + << (pass ? "pass" : std::to_string(numErrors) + " errors FAIL") + << "\n"; + return pass; +} + +/* + * Write source(ent, xi) = ent. If ownership is working, we should have the + * neighboring element with the smallest id listed on every dof holder thats + * is on a shared entity + */ +template typename Controller> +bool testOwnership( + Omega_h::Mesh &mesh, + MeshField::OmegahMeshField &omf) { + auto fieldWithCtrlr = omf.template CreateLagrangeField(); + auto field = fieldWithCtrlr.field; + const auto [shp, map] = getElement(mesh); + MeshField::FieldElement elm(mesh.nelems(), field, shp, map); + MeshField::interpolate( + elm, KOKKOS_LAMBDA(const int ent, Kokkos::Array const &) { + return Kokkos::Array{static_cast(ent)}; + }); + + MeshField::LO numErrors = 0; + // the dof holders are the vertices and, for order 2, the edges + for (int holderDim = 0; holderDim < order; ++holderDim) { + const auto topo = holderDim == 0 ? MeshField::Vertex : MeshField::Edge; + const auto up = mesh.ask_up(holderDim, dim); + const auto a2ab = up.a2ab; + const auto ab2b = up.ab2b; + MeshField::LO lerrorsDim = 0; + Kokkos::parallel_reduce( + "checkOwnership", mesh.nents(holderDim), + KOKKOS_LAMBDA(const int e, MeshField::LO &lerrors) { + auto lowest = ab2b[a2ab[e]]; + for (auto i = a2ab[e] + 1; i < a2ab[e + 1]; ++i) + lowest = Kokkos::min(lowest, ab2b[i]); + if (field(e, 0, 0, topo) != lowest) + ++lerrors; + }, + lerrorsDim); + numErrors += lerrorsDim; + } + const bool pass = numErrors == 0; + std::cout << "testOwnership(dim " << dim << ", order " << order << "): " + << (pass ? "pass" : std::to_string(numErrors) + " errors FAIL") + << "\n"; + return pass; +} + +template typename Controller> +int runTests(Omega_h::Mesh &mesh) { + MeshField::OmegahMeshField omf(mesh); + int failed = 0; + failed += !testOwnership(mesh, omf); + failed += !testOwnership(mesh, omf); + failed += !testAnalytic(mesh, omf); + failed += !testAnalytic(mesh, omf); + failed += !testAnalytic(mesh, omf); + return failed; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + int failed = 0; + { + auto lib = Omega_h::Library(&argc, &argv); + auto tris = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, + 3, 3, 0); + auto tets = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 2.0, 3.0, + 2, 2, 2); + failed += runTests<2, MeshField::KokkosController>(tris); + failed += runTests<3, MeshField::KokkosController>(tets); +#ifdef MESHFIELDS_ENABLE_CABANA + failed += runTests<2, MeshField::CabanaController>(tris); + failed += runTests<3, MeshField::CabanaController>(tets); +#endif + } + Kokkos::finalize(); + if (failed) { + std::cerr << failed << " test(s) failed\n"; + return 1; + } + std::cout << "all tests passed\n"; + return 0; +} From a4dd3117fde8bd797b4867308325827b2e6d7a33 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Tue, 6 Oct 2026 03:32:02 -0400 Subject: [PATCH 02/19] Add generic construction of coordinate fields --- src/MeshField.hpp | 83 ++++++++++++++++++++++------------ test/testIntegratorMeasure.cpp | 61 +++++++++++++++++++------ test/testInterpolate.cpp | 60 +++++++++++++++++++----- 3 files changed, 149 insertions(+), 55 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index dd541b8..9988bbb 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -32,30 +32,6 @@ MeshField::MeshInfo getMeshInfo(Omega_h::Mesh &mesh) { return meshInfo; } -template - typename Controller = MeshField::KokkosController> -decltype(MeshField::CreateCoordinateField( - MeshField::MeshInfo())) -createCoordinateField(const MeshField::MeshInfo &mesh_info, - Omega_h::Reals coords) { - const auto meshDim = mesh_info.dim; - auto coordFieldWithCtrlr = - MeshField::CreateCoordinateField( - mesh_info); - auto coordField = coordFieldWithCtrlr.field; - auto setCoordField = KOKKOS_LAMBDA(const int &i) { - coordField(i, 0, 0, MeshField::Vertex) = coords[i * meshDim]; - coordField(i, 0, 1, MeshField::Vertex) = coords[i * meshDim + 1]; - if constexpr (dim == 3) { - coordField(i, 0, 2, MeshField::Vertex) = coords[i * meshDim + 2]; - } - }; - MeshField::parallel_for(ExecutionSpace(), {0}, {mesh_info.numVtx}, - setCoordField, "setCoordField"); - return coordFieldWithCtrlr; -} - } // anonymous namespace namespace MeshField { @@ -292,15 +268,14 @@ class OmegahMeshField { Omega_h::Mesh &mesh; const MeshField::MeshInfo meshInfo; using CoordField = - decltype(createCoordinateField( - MeshField::MeshInfo(), Omega_h::Reals())); + decltype(MeshField::CreateCoordinateField(MeshField::MeshInfo())); CoordField coordField; public: OmegahMeshField(Omega_h::Mesh &mesh_in) : mesh(mesh_in), meshInfo(getMeshInfo(mesh)), - coordField(createCoordinateField( - getMeshInfo(mesh_in), mesh_in.coords())) { + coordField(CreateLagrangeCoordinateField<1>()) { static_assert(dim == 1 || dim == 2 || dim == 3); } @@ -311,8 +286,60 @@ class OmegahMeshField { order, dim, numComp>(meshInfo); } + /** + * @brief create a Lagrange coordinate field of the given order + * + * @details + * The linear field is set from the mesh coordinates. Higher order fields + * are interpolated from the linear coordinate field, so the geometry is + * straight-sided until the non-vertex nodes are moved. + * + * @tparam order the order of the coordinate field + */ + template auto CreateLagrangeCoordinateField() const { + auto fieldWithCtrlr = CreateLagrangeField(); + auto field = fieldWithCtrlr.field; + if constexpr (order == 1) { + const auto meshDim = meshInfo.dim; + const auto coords = mesh.coords(); + auto setCoordField = KOKKOS_LAMBDA(const int vtx) { + for (size_t d = 0; d < dim; ++d) + field(vtx, 0, d, MeshField::Vertex) = coords[vtx * meshDim + d]; + }; + MeshField::parallel_for(ExecutionSpace(), {0}, {meshInfo.numVtx}, + setCoordField, "setCoordField"); + } else { + const auto linear = CreateLagrangeElement(coordField.field); + auto target = CreateLagrangeElement(field); + MeshField::interpolate( + target, + KOKKOS_LAMBDA(const int ent, Kokkos::Array const &xi) { + return linear.getValue(ent, xi); + }); + } + return fieldWithCtrlr; + } + auto getCoordField() { return coordField; } + /** + * @brief create a FieldElement over the mesh elements for a Lagrange field + * + */ + template + auto CreateLagrangeElement(const Field &field) const { + static_assert(dim == 2 || dim == 3, + "CreateLagrangeElement supports triangles and tetrahedra"); + constexpr int order = Field::Order; + if constexpr (dim == 2) { + const auto [shp, map] = Omegah::getTriangleElement(mesh); + return MeshField::FieldElement(meshInfo.numTri, field, shp, map); + } else { + const auto [shp, map] = Omegah::getTetrahedronElement(mesh); + return MeshField::FieldElement(meshInfo.numTet, field, shp, map); + } + } + // FIXME support 2d and 3d and fields with order>1 template void writeVtk(Field &field) const { using FieldDataType = typename decltype(field.vtxField)::BaseType; diff --git a/test/testIntegratorMeasure.cpp b/test/testIntegratorMeasure.cpp index f2332f9..e472fe7 100644 --- a/test/testIntegratorMeasure.cpp +++ b/test/testIntegratorMeasure.cpp @@ -11,7 +11,7 @@ using ExecutionSpace = Kokkos::DefaultExecutionSpace; // compute the measure (3D volume, 2D area, 1D length) class MeasureIntegrator : public MeshField::Integrator { public: - MeasureIntegrator() : MeshField::Integrator(1) {} + MeasureIntegrator(int order) : MeshField::Integrator(order) {} void atPoints(Kokkos::View, Kokkos::View w, Kokkos::View dV) override { @@ -27,20 +27,46 @@ class MeasureIntegrator : public MeshField::Integrator { MeshField::Real measure = 0; }; -template +/* + * Move the edge nodes of the edges on y = 1 of a unit square by delta in y, + * so each of those edges becomes a parabola that adds 2/3 * length * delta + * to the area: 2/3 * delta in all. + */ +template +void moveTopEdges(Omega_h::Mesh &mesh, CoordField coordField, + MeshField::Real delta) { + const auto edgeVerts = mesh.ask_verts_of(Omega_h::EDGE); + const auto coords = mesh.coords(); + Kokkos::parallel_for( + "moveTopEdges", mesh.nedges(), KOKKOS_LAMBDA(const int edge) { + const auto a = edgeVerts[edge * 2]; + const auto b = edgeVerts[edge * 2 + 1]; + if (Kokkos::fabs(coords[a * 2 + 1] - 1) < MeshField::Epsilon && + Kokkos::fabs(coords[b * 2 + 1] - 1) < MeshField::Epsilon) + coordField(edge, 0, 1, MeshField::Edge) += delta; + }); +} + +/* + * Integrate the measure of the mesh with a coordinate field of the given + * order. + */ +template bool testMeasure(Omega_h::Mesh &mesh, MeshField::Real expected, - std::string_view name) { + std::string_view name, MeshField::Real delta = 0) { MeshField::OmegahMeshField omf(mesh); - auto coordField = omf.getCoordField().field; - const auto [shp, map] = [&]() { - if constexpr (dim == 3) { - return MeshField::Omegah::getTetrahedronElement<1>(mesh); - } else { - return MeshField::Omegah::getTriangleElement<1>(mesh); - } - }(); - MeshField::FieldElement elm(mesh.nelems(), coordField, shp, map); - MeasureIntegrator integrator; + auto coordFieldWithCtrlr = + omf.template CreateLagrangeCoordinateField(); + if constexpr (dim == 2 && order == 2) { + if (delta != 0) + moveTopEdges(mesh, coordFieldWithCtrlr.field, delta); + } else if (delta != 0) { + std::cout << "testMeasure(" << name + << "): moving edges needs a quadratic 2d field FAIL\n"; + return false; + } + auto elm = omf.CreateLagrangeElement(coordFieldWithCtrlr.field); + MeasureIntegrator integrator(order); integrator.process(elm); const auto err = std::abs(integrator.measure - expected); const bool pass = err <= MeshField::Epsilon; @@ -58,8 +84,13 @@ int main(int argc, char **argv) { 0.0, 3, 3, 0); auto tets = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 2.0, 3.0, 2, 2, 2); - failed += !testMeasure<2>(tris, 1.0, "triangles"); - failed += !testMeasure<3>(tets, 6.0, "tetrahedra"); + failed += !testMeasure<2, 1>(tris, 1.0, "triangles"); + failed += !testMeasure<3, 1>(tets, 6.0, "tetrahedra"); + failed += !testMeasure<2, 2>(tris, 1.0, "quadratic triangles"); + failed += !testMeasure<3, 2>(tets, 6.0, "quadratic tetrahedra"); + const MeshField::Real delta = 0.1; + failed += !testMeasure<2, 2>(tris, 1.0 + (2.0 / 3.0) * delta, + "curved quadratic triangles", delta); } Kokkos::finalize(); if (failed) { diff --git a/test/testInterpolate.cpp b/test/testInterpolate.cpp index ffdcb94..956b93a 100644 --- a/test/testInterpolate.cpp +++ b/test/testInterpolate.cpp @@ -35,6 +35,20 @@ KOKKOS_INLINE_FUNCTION Real linear(Kokkos::Array const &x) { return f; } +// Create a list of two parametric coordinates for evaluation +constexpr size_t numOffNodePts = 2; +template Kokkos::View offNodeCoords(size_t numElems) { + Kokkos::View lc("localCoords", numElems * numOffNodePts, dim); + Kokkos::parallel_for( + "setLocalCoords", numElems, KOKKOS_LAMBDA(const int ent) { + for (size_t d = 0; d < dim; ++d) { + lc(ent * numOffNodePts, d) = 1.0 / (dim + 1); + lc(ent * numOffNodePts + 1, d) = 0.1 * (d + 1); + } + }); + return lc; +} + /* * Interpolate an analytic function, with numComp components */ @@ -63,18 +77,9 @@ bool testAnalytic( return val; }); - // points away from the nodes: the centroid and an off-center point - constexpr size_t numPts = 2; - Kokkos::View lc("localCoords", mesh.nelems() * numPts, dim); - Kokkos::parallel_for( - "setLocalCoords", mesh.nelems(), KOKKOS_LAMBDA(const int ent) { - for (size_t d = 0; d < dim; ++d) { - lc(ent * numPts, d) = 1.0 / (dim + 1); - lc(ent * numPts + 1, d) = 0.1 * (d + 1); - } - }); - const auto values = MeshField::evaluate(elm, lc, numPts); - const auto x = MeshField::evaluate(coordElm, lc, numPts); + const auto lc = offNodeCoords(mesh.nelems()); + const auto values = MeshField::evaluate(elm, lc, numOffNodePts); + const auto x = MeshField::evaluate(coordElm, lc, numOffNodePts); MeshField::LO numErrors = 0; Kokkos::parallel_reduce( "checkAnalytic", values.extent(0), @@ -141,6 +146,36 @@ bool testOwnership( return pass; } +/* + * Ensure that straight sided quadratic coordiante field recovers + * same values as linear coordinate field + */ +template typename Controller> +bool testCoordinateField( + Omega_h::Mesh &mesh, + MeshField::OmegahMeshField &omf) { + auto quadraticWithCtrlr = omf.template CreateLagrangeCoordinateField<2>(); + auto quadratic = omf.CreateLagrangeElement(quadraticWithCtrlr.field); + auto linear = omf.CreateLagrangeElement(omf.getCoordField().field); + const auto lc = offNodeCoords(mesh.nelems()); + const auto xq = MeshField::evaluate(quadratic, lc, numOffNodePts); + const auto xl = MeshField::evaluate(linear, lc, numOffNodePts); + MeshField::LO numErrors = 0; + Kokkos::parallel_reduce( + "checkCoordinateField", xq.extent(0), + KOKKOS_LAMBDA(const int pt, MeshField::LO &lerrors) { + for (size_t d = 0; d < dim; ++d) + if (Kokkos::fabs(xq(pt, d) - xl(pt, d)) > MeshField::Epsilon) + ++lerrors; + }, + numErrors); + const bool pass = numErrors == 0; + std::cout << "testCoordinateField(dim " << dim << "): " + << (pass ? "pass" : std::to_string(numErrors) + " errors FAIL") + << "\n"; + return pass; +} + template typename Controller> int runTests(Omega_h::Mesh &mesh) { MeshField::OmegahMeshField omf(mesh); @@ -150,6 +185,7 @@ int runTests(Omega_h::Mesh &mesh) { failed += !testAnalytic(mesh, omf); failed += !testAnalytic(mesh, omf); failed += !testAnalytic(mesh, omf); + failed += !testCoordinateField(mesh, omf); return failed; } From 477e7dbd706b8c5e6cc1578cc569858f1be703f8 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 10 Oct 2026 02:10:14 -0400 Subject: [PATCH 03/19] add const to Field in MeshFieldFor --- src/MeshField_For.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/MeshField_For.hpp b/src/MeshField_For.hpp index 87a981c..3b605ad 100644 --- a/src/MeshField_For.hpp +++ b/src/MeshField_For.hpp @@ -7,7 +7,7 @@ namespace MeshField { template void parallel_for(ExecutionSpace, const std::initializer_list &start, const std::initializer_list &end, - FunctorType &vectorKernel, std::string tag) { + const FunctorType &vectorKernel, std::string tag) { constexpr auto funcRank = MeshFieldUtil::function_traits::arity; assert(funcRank >= 1); Kokkos::Array a_start = From c6d1c1c3d5856f60f1944912a980a5bdf004e92e Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 10 Oct 2026 02:10:41 -0400 Subject: [PATCH 04/19] move lambdas to functors in field creation cannot have lambdas in function with auto return type with cuda --- src/MeshField.hpp | 42 ++++++++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 9988bbb..ebd8d60 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -31,6 +31,34 @@ MeshField::MeshInfo getMeshInfo(Omega_h::Mesh &mesh) { } return meshInfo; } +// we use functors because cannot use auto return type with KOKKOS_LAMBDAS +template +struct SetCoordFieldFunctor +{ + SetCoordFieldFunctor(Field& field, Coords& coords, int meshDim) : + field_(field), coords_(coords), meshDim_(meshDim) {} + KOKKOS_INLINE_FUNCTION + void operator()(const int vtx) const { + for (size_t d = 0; d < meshDim_; ++d) + field_(vtx, 0, d, MeshField::Vertex) = coords_[vtx * meshDim_ + d]; + } + Field field_; + Coords coords_; + int meshDim_; +}; + +template +struct InterpolateElementFunctor +{ + InterpolateElementFunctor(const Element& element) : element_(element) {} + + KOKKOS_INLINE_FUNCTION + auto operator()(const int ent, Kokkos::Array const& xi) const { + return element_.getValue(ent, xi); + } + + const Element element_; +}; } // anonymous namespace @@ -299,23 +327,17 @@ class OmegahMeshField { template auto CreateLagrangeCoordinateField() const { auto fieldWithCtrlr = CreateLagrangeField(); auto field = fieldWithCtrlr.field; + // note order 1 is called in the constructor to initialize the coordinate field + // it must be called before higher order coordinate fields are constructed if constexpr (order == 1) { const auto meshDim = meshInfo.dim; const auto coords = mesh.coords(); - auto setCoordField = KOKKOS_LAMBDA(const int vtx) { - for (size_t d = 0; d < dim; ++d) - field(vtx, 0, d, MeshField::Vertex) = coords[vtx * meshDim + d]; - }; MeshField::parallel_for(ExecutionSpace(), {0}, {meshInfo.numVtx}, - setCoordField, "setCoordField"); + SetCoordFieldFunctor(field, coords, meshDim), "setCoordField"); } else { const auto linear = CreateLagrangeElement(coordField.field); auto target = CreateLagrangeElement(field); - MeshField::interpolate( - target, - KOKKOS_LAMBDA(const int ent, Kokkos::Array const &xi) { - return linear.getValue(ent, xi); - }); + MeshField::interpolate(target, InterpolateElementFunctor(linear)); } return fieldWithCtrlr; } From dc4b020893cecca8c64419c26cfb7c8d77960d18 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 10 Oct 2026 02:22:24 -0400 Subject: [PATCH 05/19] resolve review comments --- src/MeshField.hpp | 4 ++-- src/MeshField_Element.hpp | 5 +++-- test/testInterpolate.cpp | 4 ++-- 3 files changed, 7 insertions(+), 6 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index ebd8d60..176b149 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -38,7 +38,7 @@ struct SetCoordFieldFunctor SetCoordFieldFunctor(Field& field, Coords& coords, int meshDim) : field_(field), coords_(coords), meshDim_(meshDim) {} KOKKOS_INLINE_FUNCTION - void operator()(const int vtx) const { + void operator()(const size_t vtx) const { for (size_t d = 0; d < meshDim_; ++d) field_(vtx, 0, d, MeshField::Vertex) = coords_[vtx * meshDim_ + d]; } @@ -337,7 +337,7 @@ class OmegahMeshField { } else { const auto linear = CreateLagrangeElement(coordField.field); auto target = CreateLagrangeElement(field); - MeshField::interpolate(target, InterpolateElementFunctor(linear)); + MeshField::setDofByInterpolation(target, InterpolateElementFunctor(linear)); } return fieldWithCtrlr; } diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index b729b24..1ab15d7 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -651,10 +651,11 @@ DofHolderOwningEntities getDofHolderOwningEntities(FieldElement const &fes) { /** * @brief - * evaluates a source function at of the field element dof holder location + * evaluates a source function at the target field element dof holder locations + * after this function runs, the target's field values will modified */ template -void interpolate(FieldElement &target, Source const &source) { +void setDofByInterpolation(FieldElement &target, Source const &source) { constexpr size_t numNodes = decltype(target.shapeFn)::numNodes; constexpr size_t meshEntDim = FieldElement::MeshEntDim; constexpr size_t numComp = FieldElement::NumComponents; diff --git a/test/testInterpolate.cpp b/test/testInterpolate.cpp index 956b93a..8b2713e 100644 --- a/test/testInterpolate.cpp +++ b/test/testInterpolate.cpp @@ -68,7 +68,7 @@ bool testAnalytic( const auto f = KOKKOS_LAMBDA(Kokkos::Array const &x) { return order == 1 ? linear(x) : quadratic(x); }; - MeshField::interpolate( + MeshField::setDofByInterpolation( elm, KOKKOS_LAMBDA(const int ent, Kokkos::Array const &xi) { const auto x = coordElm.getValue(ent, xi); Kokkos::Array val; @@ -114,7 +114,7 @@ bool testOwnership( auto field = fieldWithCtrlr.field; const auto [shp, map] = getElement(mesh); MeshField::FieldElement elm(mesh.nelems(), field, shp, map); - MeshField::interpolate( + MeshField::setDofByInterpolation( elm, KOKKOS_LAMBDA(const int ent, Kokkos::Array const &) { return Kokkos::Array{static_cast(ent)}; }); From 4d3ca782e2f78ac2a75cd0d999add858b1193af4 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sun, 11 Oct 2026 03:53:27 -0400 Subject: [PATCH 06/19] pass functor/map into coordinate field creation --- src/MeshField.hpp | 59 ++++++++++++++++++++++++++++--------- test/testInterpolate.cpp | 63 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 108 insertions(+), 14 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 176b149..42d8201 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -32,32 +32,40 @@ MeshField::MeshInfo getMeshInfo(Omega_h::Mesh &mesh) { return meshInfo; } // we use functors because cannot use auto return type with KOKKOS_LAMBDAS -template +template struct SetCoordFieldFunctor { - SetCoordFieldFunctor(Field& field, Coords& coords, int meshDim) : - field_(field), coords_(coords), meshDim_(meshDim) {} + SetCoordFieldFunctor(Field& field, Coords& coords, int meshDim, + const CoordinateMap& map) : + field_(field), coords_(coords), meshDim_(meshDim), map_(map) {} KOKKOS_INLINE_FUNCTION void operator()(const size_t vtx) const { + Kokkos::Array x{}; for (size_t d = 0; d < meshDim_; ++d) - field_(vtx, 0, d, MeshField::Vertex) = coords_[vtx * meshDim_ + d]; + x[d] = coords_[vtx * meshDim_ + d]; + const auto mapped = map_(x); + for (size_t d = 0; d < meshDim_; ++d) + field_(vtx, 0, d, MeshField::Vertex) = mapped[d]; } Field field_; Coords coords_; int meshDim_; + const CoordinateMap map_; }; -template +template struct InterpolateElementFunctor { - InterpolateElementFunctor(const Element& element) : element_(element) {} + InterpolateElementFunctor(const Element& element, const CoordinateMap& map) + : element_(element), map_(map) {} KOKKOS_INLINE_FUNCTION auto operator()(const int ent, Kokkos::Array const& xi) const { - return element_.getValue(ent, xi); + return map_(element_.getValue(ent, xi)); } const Element element_; + const CoordinateMap map_; }; } // anonymous namespace @@ -288,6 +296,16 @@ template auto getTetrahedronElement(Omega_h::Mesh &mesh) { } // namespace Omegah +/** + * @brief the CoordinateMap that leaves the mesh geometry unchanged + */ +struct IdentityCoordinateMap { + template + KOKKOS_INLINE_FUNCTION Point operator()(Point const &x) const { + return x; + } +}; + template typename Controller = MeshField::KokkosController> @@ -318,13 +336,20 @@ class OmegahMeshField { * @brief create a Lagrange coordinate field of the given order * * @details - * The linear field is set from the mesh coordinates. Higher order fields - * are interpolated from the linear coordinate field, so the geometry is - * straight-sided until the non-vertex nodes are moved. + * Each node is placed at map(x), where x is the node's position on the + * linear (straight-sided) mesh geometry. The linear field is set from the + * mesh coordinates and higher order fields are interpolated from the + * linear coordinate field, so with the default IdentityCoordinateMap the + * geometry is the mesh's own and is straight-sided. A nonlinear map gives + * curved geometry for orders above one. * * @tparam order the order of the coordinate field + * @param map a device callable taking and returning a + * Kokkos::Array of physical coordinates */ - template auto CreateLagrangeCoordinateField() const { + template + auto CreateLagrangeCoordinateField( + const CoordinateMap &map = CoordinateMap()) const { auto fieldWithCtrlr = CreateLagrangeField(); auto field = fieldWithCtrlr.field; // note order 1 is called in the constructor to initialize the coordinate field @@ -332,12 +357,18 @@ class OmegahMeshField { if constexpr (order == 1) { const auto meshDim = meshInfo.dim; const auto coords = mesh.coords(); - MeshField::parallel_for(ExecutionSpace(), {0}, {meshInfo.numVtx}, - SetCoordFieldFunctor(field, coords, meshDim), "setCoordField"); + MeshField::parallel_for( + ExecutionSpace(), {0}, {meshInfo.numVtx}, + SetCoordFieldFunctor(field, coords, meshDim, map), + "setCoordField"); } else { const auto linear = CreateLagrangeElement(coordField.field); auto target = CreateLagrangeElement(field); - MeshField::setDofByInterpolation(target, InterpolateElementFunctor(linear)); + MeshField::setDofByInterpolation( + target, + InterpolateElementFunctor( + linear, map)); } return fieldWithCtrlr; } diff --git a/test/testInterpolate.cpp b/test/testInterpolate.cpp index 8b2713e..7d365d0 100644 --- a/test/testInterpolate.cpp +++ b/test/testInterpolate.cpp @@ -176,6 +176,66 @@ bool testCoordinateField( return pass; } +// affine, so coordinate fields of any order reproduce it exactly +template struct Shear { + KOKKOS_INLINE_FUNCTION Kokkos::Array + operator()(Kokkos::Array x) const { + x[0] += 0.3 * x[1]; + if constexpr (dim == 3) { + x[0] += 0.1 * x[2]; + x[1] += 0.2 * x[2]; + } + return x; + } +}; + +// quadratic, so quadratic coordinate fields reproduce it exactly +template struct Bend { + KOKKOS_INLINE_FUNCTION Kokkos::Array + operator()(Kokkos::Array x) const { + x[1] += 0.2 * x[0] * x[0]; + return x; + } +}; + +/* + * Ensure that a coordinate field created with a map recovers the map applied + * to the linear coordinate field + */ +template typename Controller> +bool testCoordinateMap( + Omega_h::Mesh &mesh, + MeshField::OmegahMeshField &omf, + CoordinateMap const &coordMap, std::string const &name) { + auto mappedWithCtrlr = + omf.template CreateLagrangeCoordinateField(coordMap); + auto mapped = omf.CreateLagrangeElement(mappedWithCtrlr.field); + auto linear = omf.CreateLagrangeElement(omf.getCoordField().field); + const auto lc = offNodeCoords(mesh.nelems()); + const auto xm = MeshField::evaluate(mapped, lc, numOffNodePts); + const auto xl = MeshField::evaluate(linear, lc, numOffNodePts); + MeshField::LO numErrors = 0; + Kokkos::parallel_reduce( + "checkCoordinateMap", xm.extent(0), + KOKKOS_LAMBDA(const int pt, MeshField::LO &lerrors) { + Kokkos::Array x; + for (size_t d = 0; d < dim; ++d) + x[d] = xl(pt, d); + const auto expected = coordMap(x); + for (size_t d = 0; d < dim; ++d) + if (Kokkos::fabs(xm(pt, d) - expected[d]) > MeshField::Epsilon) + ++lerrors; + }, + numErrors); + const bool pass = numErrors == 0; + std::cout << "testCoordinateMap(dim " << dim << ", order " << order << ", " + << name << "): " + << (pass ? "pass" : std::to_string(numErrors) + " errors FAIL") + << "\n"; + return pass; +} + template typename Controller> int runTests(Omega_h::Mesh &mesh) { MeshField::OmegahMeshField omf(mesh); @@ -186,6 +246,9 @@ int runTests(Omega_h::Mesh &mesh) { failed += !testAnalytic(mesh, omf); failed += !testAnalytic(mesh, omf); failed += !testCoordinateField(mesh, omf); + failed += !testCoordinateMap(mesh, omf, Shear(), "shear"); + failed += !testCoordinateMap(mesh, omf, Shear(), "shear"); + failed += !testCoordinateMap(mesh, omf, Bend(), "bend"); return failed; } From bec167865de24f2a3b3e1af7e21ebdddd0c912f8 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sun, 11 Oct 2026 13:19:03 -0400 Subject: [PATCH 07/19] add shared test functions header --- test/MeshField_TestFunctions.hpp | 83 ++++++++++++++++++++++++++++++++ test/testIntegratorMeasure.cpp | 52 +++++++------------- test/testInterpolate.cpp | 41 +--------------- test/testOmegahTet.cpp | 41 +++++----------- test/testOmegahTri.cpp | 39 +++++---------- 5 files changed, 129 insertions(+), 127 deletions(-) create mode 100644 test/MeshField_TestFunctions.hpp diff --git a/test/MeshField_TestFunctions.hpp b/test/MeshField_TestFunctions.hpp new file mode 100644 index 0000000..e0997a3 --- /dev/null +++ b/test/MeshField_TestFunctions.hpp @@ -0,0 +1,83 @@ +#ifndef MESHFIELD_TESTFUNCTIONS_HPP +#define MESHFIELD_TESTFUNCTIONS_HPP + +/** + * @file MeshField_TestFunctions.hpp + * @brief analytic functions and coordinate maps shared by the tests + */ + +#include +#include + +template using Vec = Kokkos::Array; + +/** @brief f(x) = 1 + sum_d (d+1) x_d */ +template struct LinearFunction { + KOKKOS_INLINE_FUNCTION MeshField::Real operator()(Vec const &x) const { + MeshField::Real f = 1; + for (size_t d = 0; d < dim; ++d) + f += (d + 1) * x[d]; + return f; + } + KOKKOS_INLINE_FUNCTION Vec grad(Vec const &) const { + Vec g; + for (size_t d = 0; d < dim; ++d) + g[d] = d + 1; + return g; + } +}; + +/** + * @brief f(x) = 1 + sum_d ((d+1) x_d + x_d^2) + x_0 x_1 + * + * The cross term gives the Hessian an off-diagonal entry. + */ +template struct QuadraticFunction { + static_assert(dim >= 2, "the cross term needs at least two dimensions"); + KOKKOS_INLINE_FUNCTION MeshField::Real operator()(Vec const &x) const { + MeshField::Real f = 1 + x[0] * x[1]; + for (size_t d = 0; d < dim; ++d) + f += (d + 1) * x[d] + x[d] * x[d]; + return f; + } + KOKKOS_INLINE_FUNCTION Vec grad(Vec const &x) const { + Vec g; + for (size_t d = 0; d < dim; ++d) + g[d] = (d + 1) + 2 * x[d]; + g[0] += x[1]; + g[1] += x[0]; + return g; + } +}; + +/** + * @brief CoordinateMap for an affine shear that makes every Jacobian a full + * matrix: x += 0.3y + 0.1z, y += 0.2z + * + * Its determinant is one, so it preserves measure. Coordinate fields of any + * order reproduce it exactly. + */ +template struct Shear { + KOKKOS_INLINE_FUNCTION Vec operator()(Vec x) const { + x[0] += 0.3 * x[1]; + if constexpr (dim == 3) { + x[0] += 0.1 * x[2]; + x[1] += 0.2 * x[2]; + } + return x; + } +}; + +/** + * @brief CoordinateMap that bends the mesh: y += 0.2x^2 + * + * Quadratic, so quadratic coordinate fields reproduce it exactly. + */ +template struct Bend { + KOKKOS_INLINE_FUNCTION Vec operator()(Vec x) const { + x[1] += 0.2 * x[0] * x[0]; + return x; + } +}; + +#endif diff --git a/test/testIntegratorMeasure.cpp b/test/testIntegratorMeasure.cpp index e472fe7..04a4f6d 100644 --- a/test/testIntegratorMeasure.cpp +++ b/test/testIntegratorMeasure.cpp @@ -2,6 +2,7 @@ #include "MeshField.hpp" #include "MeshField_Element.hpp" #include "MeshField_Integrate.hpp" +#include "MeshField_TestFunctions.hpp" #include "Omega_h_build.hpp" #include #include @@ -27,44 +28,18 @@ class MeasureIntegrator : public MeshField::Integrator { MeshField::Real measure = 0; }; -/* - * Move the edge nodes of the edges on y = 1 of a unit square by delta in y, - * so each of those edges becomes a parabola that adds 2/3 * length * delta - * to the area: 2/3 * delta in all. - */ -template -void moveTopEdges(Omega_h::Mesh &mesh, CoordField coordField, - MeshField::Real delta) { - const auto edgeVerts = mesh.ask_verts_of(Omega_h::EDGE); - const auto coords = mesh.coords(); - Kokkos::parallel_for( - "moveTopEdges", mesh.nedges(), KOKKOS_LAMBDA(const int edge) { - const auto a = edgeVerts[edge * 2]; - const auto b = edgeVerts[edge * 2 + 1]; - if (Kokkos::fabs(coords[a * 2 + 1] - 1) < MeshField::Epsilon && - Kokkos::fabs(coords[b * 2 + 1] - 1) < MeshField::Epsilon) - coordField(edge, 0, 1, MeshField::Edge) += delta; - }); -} - /* * Integrate the measure of the mesh with a coordinate field of the given - * order. + * order, created with coordMap. */ -template +template bool testMeasure(Omega_h::Mesh &mesh, MeshField::Real expected, - std::string_view name, MeshField::Real delta = 0) { + std::string_view name, + CoordinateMap const &coordMap = CoordinateMap()) { MeshField::OmegahMeshField omf(mesh); auto coordFieldWithCtrlr = - omf.template CreateLagrangeCoordinateField(); - if constexpr (dim == 2 && order == 2) { - if (delta != 0) - moveTopEdges(mesh, coordFieldWithCtrlr.field, delta); - } else if (delta != 0) { - std::cout << "testMeasure(" << name - << "): moving edges needs a quadratic 2d field FAIL\n"; - return false; - } + omf.template CreateLagrangeCoordinateField(coordMap); auto elm = omf.CreateLagrangeElement(coordFieldWithCtrlr.field); MeasureIntegrator integrator(order); integrator.process(elm); @@ -88,9 +63,18 @@ int main(int argc, char **argv) { failed += !testMeasure<3, 1>(tets, 6.0, "tetrahedra"); failed += !testMeasure<2, 2>(tris, 1.0, "quadratic triangles"); failed += !testMeasure<3, 2>(tets, 6.0, "quadratic tetrahedra"); + failed += !testMeasure<3, 2>(tets, 6.0, "sheared quadratic tetrahedra", + Shear<3>()); + // the unit square with its top bent into y = 1 + delta * x(1 - x). The + // top edges are parabolas in their edge parameter, so quadratic elements + // reproduce the boundary, and the area, exactly. const MeshField::Real delta = 0.1; - failed += !testMeasure<2, 2>(tris, 1.0 + (2.0 / 3.0) * delta, - "curved quadratic triangles", delta); + failed += !testMeasure<2, 2>( + tris, 1.0 + delta / 6.0, "curved quadratic triangles", + KOKKOS_LAMBDA(Vec<2> x) { + x[1] *= 1 + delta * x[0] * (1 - x[0]); + return x; + }); } Kokkos::finalize(); if (failed) { diff --git a/test/testInterpolate.cpp b/test/testInterpolate.cpp index 7d365d0..71233f1 100644 --- a/test/testInterpolate.cpp +++ b/test/testInterpolate.cpp @@ -1,5 +1,6 @@ #include "KokkosController.hpp" #include "MeshField.hpp" +#include "MeshField_TestFunctions.hpp" #ifdef MESHFIELDS_ENABLE_CABANA #include "CabanaController.hpp" #endif @@ -19,22 +20,6 @@ template auto getElement(Omega_h::Mesh &mesh) { } } -template -KOKKOS_INLINE_FUNCTION Real quadratic(Kokkos::Array const &x) { - Real f = 1; - for (size_t d = 0; d < dim; ++d) - f += (d + 1) * x[d] + x[d] * x[d]; - return f; -} - -template -KOKKOS_INLINE_FUNCTION Real linear(Kokkos::Array const &x) { - Real f = 1; - for (size_t d = 0; d < dim; ++d) - f += (d + 1) * x[d]; - return f; -} - // Create a list of two parametric coordinates for evaluation constexpr size_t numOffNodePts = 2; template Kokkos::View offNodeCoords(size_t numElems) { @@ -66,7 +51,7 @@ bool testAnalytic( MeshField::FieldElement elm(mesh.nelems(), fieldWithCtrlr.field, shp, map); const auto f = KOKKOS_LAMBDA(Kokkos::Array const &x) { - return order == 1 ? linear(x) : quadratic(x); + return order == 1 ? LinearFunction()(x) : QuadraticFunction()(x); }; MeshField::setDofByInterpolation( elm, KOKKOS_LAMBDA(const int ent, Kokkos::Array const &xi) { @@ -176,28 +161,6 @@ bool testCoordinateField( return pass; } -// affine, so coordinate fields of any order reproduce it exactly -template struct Shear { - KOKKOS_INLINE_FUNCTION Kokkos::Array - operator()(Kokkos::Array x) const { - x[0] += 0.3 * x[1]; - if constexpr (dim == 3) { - x[0] += 0.1 * x[2]; - x[1] += 0.2 * x[2]; - } - return x; - } -}; - -// quadratic, so quadratic coordinate fields reproduce it exactly -template struct Bend { - KOKKOS_INLINE_FUNCTION Kokkos::Array - operator()(Kokkos::Array x) const { - x[1] += 0.2 * x[0] * x[0]; - return x; - } -}; - /* * Ensure that a coordinate field created with a map recovers the map applied * to the linear coordinate field diff --git a/test/testOmegahTet.cpp b/test/testOmegahTet.cpp index bb2a500..2514e2b 100644 --- a/test/testOmegahTet.cpp +++ b/test/testOmegahTet.cpp @@ -1,5 +1,6 @@ #include "KokkosController.hpp" #include "MeshField.hpp" +#include "MeshField_TestFunctions.hpp" #include "MeshField_Element.hpp" //remove? #include "MeshField_Fail.hpp" //remove? #include "MeshField_For.hpp" //remove? @@ -17,22 +18,6 @@ using ExecutionSpace = Kokkos::DefaultExecutionSpace; using MemorySpace = Kokkos::DefaultExecutionSpace::memory_space; -struct LinearFunction { - KOKKOS_INLINE_FUNCTION - MeshField::Real operator()(MeshField::Real x, MeshField::Real y, - MeshField::Real z) const { - return 2.0 * x + y + z; - } -}; - -struct QuadraticFunction { - KOKKOS_INLINE_FUNCTION - MeshField::Real operator()(MeshField::Real x, MeshField::Real y, - MeshField::Real z) const { - return (x * x) + (2.0 * y) + z; - } -}; - Omega_h::Mesh createMeshTet(Omega_h::Library &lib) { auto world = lib.world(); const auto family = OMEGA_H_SIMPLEX; @@ -66,12 +51,12 @@ bool checkResult(Omega_h::Mesh &mesh, Result &result, CoordField coordField, const auto x = globalCoords(pt, 0); const auto y = globalCoords(pt, 1); const auto z = globalCoords(pt, 2); - const auto expected = func(x, y, z); + const auto expected = func(Vec<3>{x, y, z}); for (size_t i = 0; i < numComp; ++i) { const auto computed = result(pt, i); MeshField::LO isError = 0; if (Kokkos::fabs(computed - expected) > - MeshField::MachinePrecision) { + MeshField::Epsilon) { isError = 1; Kokkos::printf( "result for elm %d, pt %d, does not match: expected " @@ -97,7 +82,7 @@ void setVertices(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) { const auto y = coords[vtx * MeshDim + 1]; const auto z = coords[vtx * MeshDim + 2]; for (size_t i = 0; i < field.numComp; ++i) { - field(vtx, 0, i, MeshField::Vertex) = func(x, y, z); + field(vtx, 0, i, MeshField::Vertex) = func(Vec<3>{x, y, z}); } }; MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()}, @@ -122,7 +107,7 @@ void setEdges(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) { const auto z = (coords[left * MeshDim + 2] + coords[right * MeshDim + 2]) / 2.0; for (size_t i = 0; i < field.numComp; ++i) { - field(edge, 0, i, MeshField::Edge) = func(x, y, z); + field(edge, 0, i, MeshField::Edge) = func(Vec<3>{x, y, z}); } }; MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nedges()}, @@ -180,7 +165,7 @@ void runTest(Omega_h::Mesh &mesh, if (failed) { std::string fieldErr = ShapeOrder == 1 ? "linear" : "quadratic"; std::string functionErr; - if constexpr (std::is_same_v) { + if constexpr (std::is_same_v>) { functionErr = "linear"; } else { functionErr = "quadratic"; @@ -237,42 +222,42 @@ void doRun(Omega_h::Mesh &mesh, { const auto ShapeOrder = 1; const auto numComponents = 1; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<3>()); } { const auto ShapeOrder = 2; const auto numComponents = 1; runTest(mesh, omf, testCase, - QuadraticFunction()); + QuadraticFunction<3>()); } { const auto ShapeOrder = 2; const auto numComponents = 1; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<3>()); } { const auto ShapeOrder = 1; const auto numComponents = 2; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<3>()); } { const auto ShapeOrder = 1; const auto numComponents = 3; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<3>()); } { const auto ShapeOrder = 2; const auto numComponents = 2; runTest(mesh, omf, testCase, - QuadraticFunction()); + QuadraticFunction<3>()); } { const auto ShapeOrder = 2; const auto numComponents = 3; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<3>()); } } } diff --git a/test/testOmegahTri.cpp b/test/testOmegahTri.cpp index 09fc68e..c965de1 100644 --- a/test/testOmegahTri.cpp +++ b/test/testOmegahTri.cpp @@ -1,5 +1,6 @@ #include "KokkosController.hpp" #include "MeshField.hpp" +#include "MeshField_TestFunctions.hpp" #include "MeshField_Element.hpp" //remove? #include "MeshField_Fail.hpp" //remove? #include "MeshField_For.hpp" //remove? @@ -17,20 +18,6 @@ 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; - } -}; - -struct QuadraticFunction { - KOKKOS_INLINE_FUNCTION - MeshField::Real operator()(MeshField::Real x, MeshField::Real y) const { - return (x * x) + (2.0 * y); - } -}; - Omega_h::Mesh createMeshTri18(Omega_h::Library &lib) { auto world = lib.world(); const auto family = OMEGA_H_SIMPLEX; @@ -63,12 +50,12 @@ bool checkResult(Omega_h::Mesh &mesh, Result &result, CoordField coordField, 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); + const auto expected = func(Vec<2>{x, y}); for (size_t i = 0; i < numComp; ++i) { const auto computed = result(pt, i); MeshField::LO isError = 0; if (Kokkos::fabs(computed - expected) > - MeshField::MachinePrecision) { + MeshField::Epsilon) { isError = 1; Kokkos::printf( "result for elm %d, pt %d, does not match: expected " @@ -93,7 +80,7 @@ void setVertices(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) { const auto x = coords[vtx * MeshDim]; const auto y = coords[vtx * MeshDim + 1]; for (size_t i = 0; i < field.numComp; ++i) { - field(vtx, 0, i, MeshField::Vertex) = func(x, y); + field(vtx, 0, i, MeshField::Vertex) = func(Vec<2>{x, y}); } }; MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nverts()}, @@ -116,7 +103,7 @@ void setEdges(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) { const auto y = (coords[left * MeshDim + 1] + coords[right * MeshDim + 1]) / 2.0; for (size_t i = 0; i < field.numComp; ++i) { - field(edge, 0, i, MeshField::Edge) = func(x, y); + field(edge, 0, i, MeshField::Edge) = func(Vec<2>{x, y}); } }; MeshField::parallel_for(ExecutionSpace(), {0}, {mesh.nedges()}, @@ -173,7 +160,7 @@ void runTest(Omega_h::Mesh &mesh, if (failed) { std::string fieldErr = ShapeOrder == 1 ? "linear" : "quadratic"; std::string functionErr; - if constexpr (std::is_same_v) { + if constexpr (std::is_same_v>) { functionErr = "linear"; } else { functionErr = "quadratic"; @@ -223,43 +210,43 @@ void doRun(Omega_h::Mesh &mesh, { const auto ShapeOrder = 1; const auto numComponents = 1; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<2>()); } { const auto ShapeOrder = 2; const auto numComponents = 1; runTest(mesh, omf, testCase, - QuadraticFunction()); + QuadraticFunction<2>()); } { const auto ShapeOrder = 2; const auto numComponents = 1; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<2>()); } { const auto ShapeOrder = 1; const auto numComponents = 2; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<2>()); } { const auto ShapeOrder = 1; const auto numComponents = 3; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<2>()); } { const auto ShapeOrder = 2; const auto numComponents = 2; runTest(mesh, omf, testCase, - QuadraticFunction()); + QuadraticFunction<2>()); } { const auto ShapeOrder = 2; const auto numComponents = 3; - runTest(mesh, omf, testCase, LinearFunction()); + runTest(mesh, omf, testCase, LinearFunction<2>()); } } } From e20050118e66c6f814434a110c05e2c1b89a0fc1 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sun, 11 Oct 2026 02:25:58 -0400 Subject: [PATCH 08/19] split coefficients and value computation To support non-isoparametic elements, we split the evaluation of the nodal coefficients and the multiplication with the shape functions --- src/MeshField_Element.hpp | 132 ++++++++++++++++++++++++++++++++------ 1 file changed, 111 insertions(+), 21 deletions(-) diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index 1ab15d7..93d8b6e 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -173,40 +173,129 @@ struct FieldElement { FieldAccessor::numComp>; static const size_t NumComponents = FieldAccessor::numComp; + /** @brief every dof of one element, indexed [node * NumComponents + comp] */ + using DofArray = + Kokkos::Array::type, + ShapeType::numNodes * NumComponents>; + + /** @brief parametric gradients of every component at one point, indexed + * [direction * NumComponents + comp] */ + using LocalGradArray = + Kokkos::Array::type, + MeshEntDim * NumComponents>; + /** * @brief - * evaluate the field in the specified element at the specified - * parametric/local/area coordinate + * gather element's dofs into an array * - * @todo add expression in documetation for summation of shape function * - * field values - * - * @details - * heavily based on SCOREC/core @ 7cd76473 apf/apfElement.cc * * @param ent the mesh entity index - * @param localCoord the parametric coordinate - * @return the result of evaluation + * @return the element's dofs, indexed [node * NumComponents + component] */ - KOKKOS_INLINE_FUNCTION ValArray - getValue(int ent, Kokkos::Array localCoord) const { + KOKKOS_INLINE_FUNCTION DofArray gatherDofs(int ent) const { assert(ent >= 0); assert(static_cast(ent) < numMeshEnts); - ValArray c; - const auto shapeValues = shapeFn.getValues(localCoord); - for (size_t ci = 0; ci < NumComponents; ++ci) - c[ci] = 0; + DofArray dofs; for (auto topo : elm2dof.getTopology()) { // element topology - for (size_t ni = 0; ni < shapeFn.numNodes; ++ni) { + for (size_t ni = 0; ni < ShapeType::numNodes; ++ni) { for (size_t ci = 0; ci < NumComponents; ++ci) { - auto map = elm2dof(ni, ci, ent, topo); - const auto fval = + const auto map = elm2dof(ni, ci, ent, topo); + dofs[ni * NumComponents + ci] = field(map.entity, map.node, map.component, map.topo); - c[ci] += fval * shapeValues[ni]; } } } - return c; + return dofs; + } + + /** + * @brief + * the coefficients of one element on its reference basis + * + * @details + * The coefficients depend only on the element, not on the evaluation point, + * so compute them once per element and pass them to valueFromCoefficients + * and localGradientFromCoefficients for each point. For a FieldElement the + * stored dofs are the coefficients; elements whose dofs must be transformed + * to the reference element (e.g., the Bell element) override this. + * + * @param ent the mesh entity index + * @return coefficients, indexed [node * NumComponents + component] + */ + KOKKOS_INLINE_FUNCTION DofArray coefficients(int ent) const { + return gatherDofs(ent); + } + + /** + * @brief + * contract an element's coefficients with the basis at a parametric + * coordinate + * + * @details + * heavily based on SCOREC/core @ 7cd76473 apf/apfElement.cc + * + * @param c coefficients from coefficients(ent) + * @param localCoord the parametric coordinate + * @return the field value + */ + template + KOKKOS_INLINE_FUNCTION ValArray valueFromCoefficients( + Coefficients const &c, Kokkos::Array localCoord) const { + const auto shapeValues = shapeFn.getValues(localCoord); + ValArray v; + for (size_t ci = 0; ci < NumComponents; ++ci) + v[ci] = 0; + for (size_t ni = 0; ni < ShapeType::numNodes; ++ni) + for (size_t ci = 0; ci < NumComponents; ++ci) + v[ci] += c[ni * NumComponents + ci] * shapeValues[ni]; + return v; + } + + /** + * @brief + * contract an element's coefficients with the basis gradients at a + * parametric coordinate + * + * @details + * The result is the gradient with respect to the parametric coordinates, + * not the physical coordinates; mapping it to the physical element requires + * the geometry's Jacobian. + * + * @param c coefficients from coefficients(ent) + * @param localCoord the parametric coordinate + * @return parametric gradients, indexed [direction * NumComponents + component] + */ + template + KOKKOS_INLINE_FUNCTION LocalGradArray localGradientFromCoefficients( + Coefficients const &c, Kokkos::Array localCoord) const { + const auto shapeGrads = shapeFn.getLocalGradients(localCoord); + LocalGradArray g; + for (size_t i = 0; i < MeshEntDim * NumComponents; ++i) + g[i] = 0; + for (size_t ni = 0; ni < ShapeType::numNodes; ++ni) + for (size_t d = 0; d < MeshEntDim; ++d) + for (size_t ci = 0; ci < NumComponents; ++ci) + g[d * NumComponents + ci] += + c[ni * NumComponents + ci] * shapeGrads[ni * MeshEntDim + d]; + return g; + } + + /** + * @brief + * evaluate the field in the specified element at the specified + * parametric/local/area coordinate + * + * @details + * Prefer coefficients and valueFromCoefficients when evaluating an element + * at more than one point. + * + * @param ent the mesh entity index + * @param localCoord the parametric coordinate + * @return the result of evaluation + */ + KOKKOS_INLINE_FUNCTION ValArray + getValue(int ent, Kokkos::Array localCoord) const { + return valueFromCoefficients(coefficients(ent), localCoord); } using NodeArray = @@ -516,11 +605,12 @@ evaluate(FieldElement &fes, Kokkos::View localCoords, Kokkos::parallel_for( fes.numMeshEnts, KOKKOS_LAMBDA(const int ent) { Kokkos::Array lc; + const auto c = fes.coefficients(ent); // TODO use nested parallel for? for (auto pt = offsets(ent); pt < offsets(ent + 1); pt++) { for (size_t i = 0; i < localCoords.extent(1); i++) // better way? lc[i] = localCoords(pt, i); - const auto val = fes.getValue(ent, lc); + const auto val = fes.valueFromCoefficients(c, lc); for (size_t i = 0; i < numComponents; i++) res(pt, i) = val[i]; } From caa297b08ac5cdfe2634d64a5e0b20abc476c274 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 26 Sep 2026 18:30:26 -0400 Subject: [PATCH 09/19] add MappedElement for fields on a separate geometry element MappedElement composes a field element with a geometry element over the coordinate field (see #100) --- CMakeLists.txt | 4 + src/MeshField.hpp | 18 +++ src/MeshField_Element.hpp | 37 +++++ src/MeshField_Integrate.hpp | 2 +- src/MeshField_LinearAlgebra.hpp | 72 ++++++++++ src/MeshField_MappedElement.hpp | 197 ++++++++++++++++++++++++++ test/testOmegahMappedElement.cpp | 231 +++++++++++++++++++++++++++++++ 7 files changed, 560 insertions(+), 1 deletion(-) create mode 100644 src/MeshField_LinearAlgebra.hpp create mode 100644 src/MeshField_MappedElement.hpp create mode 100644 test/testOmegahMappedElement.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index f295708..2fd18bb 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -26,6 +26,8 @@ set(MESHFIELD_HEADERS src/MeshField_Macros.hpp src/KokkosController.hpp src/MeshField_Element.hpp + src/MeshField_LinearAlgebra.hpp + src/MeshField_MappedElement.hpp src/MeshField_Integrate.hpp src/MeshField_Shape.hpp src/MeshField_ShapeField.hpp @@ -214,6 +216,7 @@ meshfields_add_exe(ExceptionTest test/testExceptions.cpp) meshfields_add_exe(OmegahTetTest test/testOmegahTet.cpp) meshfields_add_exe(IntegratorMeasure test/testIntegratorMeasure.cpp) meshfields_add_exe(Interpolate test/testInterpolate.cpp) +meshfields_add_exe(OmegahMappedElement test/testOmegahMappedElement.cpp) meshfields_add_exe(OmegahThwaitesSprAdapt test/testSprThwaitesAdapt.cpp) if(MeshFields_USE_Cabana) @@ -237,6 +240,7 @@ smoke_test_func(OmegahTriTests ./OmegahTriTests) test_func(OmegahTetTest ./OmegahTetTest) test_func(IntegratorMeasure ./IntegratorMeasure) test_func(Interpolate ./Interpolate) +test_func(OmegahMappedElement ./OmegahMappedElement) if(MeshFields_USE_EXCEPTIONS) # exception caught - no error test_func(ExceptionTest ./ExceptionTest) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 42d8201..5202ff2 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -5,6 +5,7 @@ #include "MeshField_Element.hpp" #include "MeshField_Fail.hpp" #include "MeshField_For.hpp" +#include "MeshField_MappedElement.hpp" #include "MeshField_ShapeField.hpp" #include "Omega_h_file.hpp" //move #include "Omega_h_mesh.hpp" //move @@ -378,6 +379,10 @@ class OmegahMeshField { /** * @brief create a FieldElement over the mesh elements for a Lagrange field * + * @details + * The shape and dof map are chosen from the field's order. + * + * @param coordinate field used to create the element */ template auto CreateLagrangeElement(const Field &field) const { @@ -393,6 +398,19 @@ class OmegahMeshField { } } + /** + * @brief create a MappedElement for a field on the given geometry + * + * @param field the Lagrange field to evaluate + * @param geomElem the geometry element: CreateLagrangeElement of a + * coordinate field + */ + template + auto CreateLagrangeElement(const ShapeField &field, + const GeomElem &geomElem) const { + return MeshField::MappedElement(CreateLagrangeElement(field), geomElem); + } + // FIXME support 2d and 3d and fields with order>1 template void writeVtk(Field &field) const { using FieldDataType = typename decltype(field.vtxField)::BaseType; diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index 93d8b6e..341b7b4 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -157,6 +157,12 @@ struct FieldElement { const ElementDofHolderAccessor elm2dofIn) : numMeshEnts(numMeshEntsIn), field(fieldIn), shapeFn(shapeFnIn), elm2dof(elm2dofIn) {} + + /** @brief the topology of the element's reference element */ + static constexpr KOKKOS_FUNCTION auto getTopology() { + return ElementDofHolderAccessor::getTopology(); + } + /* general template for baseType which simply sets type */ template struct baseType { @@ -316,6 +322,37 @@ struct FieldElement { return c; } + /** @brief the Jacobian at one point, row-major with J[i][j] = + * dx_j/dxi_i */ + using JacobianArray = Kokkos::Array; + + /** + * @brief + * compute the Jacobian at a parametric coordinate from the node coordinates + * + * @details + * The field must be a coordinate field. + * + * heavily based on SCOREC/core @ 7cd76473 apf/apfVectorElement.cc + * + * @param x node coordinates from getNodeValues + * @param localCoord the parametric coordinate + * @return the Jacobian, row-major with J[i][j] = dx_j/dxi_i + */ + KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromNodeValues( + NodeArray const &x, Kokkos::Array localCoord) const { + const auto grads = shapeFn.getLocalGradients(localCoord); + JacobianArray J; + for (size_t i = 0; i < MeshEntDim * MeshEntDim; ++i) + J[i] = 0; + for (size_t n = 0; n < ShapeType::numNodes; ++n) + for (size_t i = 0; i < MeshEntDim; ++i) + for (size_t j = 0; j < MeshEntDim; ++j) + J[i * MeshEntDim + j] += + grads[n * MeshEntDim + i] * x[n * MeshEntDim + j]; + return J; + } + /** * @brief * compute the Jacobian of an edge diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index da80ba4..04bd759 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -253,7 +253,7 @@ class Integrator { * \param fes FieldElement * */ template void process(FieldElement &fes) { - constexpr auto topo = decltype(FieldElement::elm2dof)::getTopology(); + constexpr auto topo = FieldElement::getTopology(); pre(); auto ip = getIntegrationPoints(order); auto localCoords = getIntegrationPointLocalCoords(fes, ip); diff --git a/src/MeshField_LinearAlgebra.hpp b/src/MeshField_LinearAlgebra.hpp new file mode 100644 index 0000000..790c10d --- /dev/null +++ b/src/MeshField_LinearAlgebra.hpp @@ -0,0 +1,72 @@ +#ifndef MESHFIELD_LINEARALGEBRA_HPP +#define MESHFIELD_LINEARALGEBRA_HPP + +#include +#include +#include + +/** + * @file MeshField_LinearAlgebra.hpp + * @brief small dense matrix operations callable on host and device + * + * Matrices are Dim x Dim and stored row-major in a + * Kokkos::Array. + */ + +namespace MeshField { + +/** + * @brief the determinant of a 1x1, 2x2, or 3x3 matrix stored row-major + */ +template +KOKKOS_INLINE_FUNCTION Real +determinant(Kokkos::Array const &A) { + static_assert(Dim >= 1 && Dim <= 3, "only 1x1, 2x2, and 3x3 are supported"); + if constexpr (Dim == 1) { + return A[0]; + } else if constexpr (Dim == 2) { + return A[0] * A[3] - A[1] * A[2]; + } else { + return A[0] * (A[4] * A[8] - A[5] * A[7]) - + A[1] * (A[3] * A[8] - A[5] * A[6]) + + A[2] * (A[3] * A[7] - A[4] * A[6]); + } +} + +/** + * @brief invert a 1x1, 2x2, or 3x3 matrix stored row-major + */ +template +KOKKOS_INLINE_FUNCTION Kokkos::Array +invert(Kokkos::Array const &A) { + static_assert(Dim >= 1 && Dim <= 3, "only 1x1, 2x2, and 3x3 are supported"); + const auto det = determinant(A); + assert(Kokkos::fabs(det) > Epsilon); + Kokkos::Array inv; + if constexpr (Dim == 1) { + inv[0] = 1; + } else if constexpr (Dim == 2) { + inv[0] = A[3]; + inv[1] = -A[1]; + inv[2] = -A[2]; + inv[3] = A[0]; + } else { + // adjugate + inv[0] = A[4] * A[8] - A[5] * A[7]; + inv[1] = A[2] * A[7] - A[1] * A[8]; + inv[2] = A[1] * A[5] - A[2] * A[4]; + inv[3] = A[5] * A[6] - A[3] * A[8]; + inv[4] = A[0] * A[8] - A[2] * A[6]; + inv[5] = A[2] * A[3] - A[0] * A[5]; + inv[6] = A[3] * A[7] - A[4] * A[6]; + inv[7] = A[1] * A[6] - A[0] * A[7]; + inv[8] = A[0] * A[4] - A[1] * A[3]; + } + for (size_t i = 0; i < Dim * Dim; ++i) + inv[i] /= det; + return inv; +} + +} // namespace MeshField + +#endif diff --git a/src/MeshField_MappedElement.hpp b/src/MeshField_MappedElement.hpp new file mode 100644 index 0000000..8c5477e --- /dev/null +++ b/src/MeshField_MappedElement.hpp @@ -0,0 +1,197 @@ +#ifndef MESHFIELD_MAPPEDELEMENT_HPP +#define MESHFIELD_MAPPEDELEMENT_HPP + +#include +#include +#include +#include +#include +#include // std::declval + +/** + * @file MeshField_MappedElement.hpp + * @brief Elements whose geometry, dofs, or values need mapping between the + * reference element and the physical element. + * + * A MappedElement composes three things: + * + * - **a geometry element**: a FieldElement over the coordinate field that + * supplies the map from the reference element to the physical element. + * - **a DofTransformation**: converts the dofs stored on mesh entities to + * coefficients on the reference basis. See: + * R.C. Kirby, "A general approach to transforming finite + * elements", SMAI J. Comput. Math. 4 (2018) 197-224 , and https://defelement.org/ + * - **a MapType**: converts values and derivatives on the reference element + * to the physical element. See: https://defelement.org + */ + +namespace MeshField { + +struct IdentityDofTransformation { + template + KOKKOS_INLINE_FUNCTION auto coefficients(FieldElem const &field, + GeomElem const &, int ent) const { + return field.coefficients(ent); + } +}; + +struct IdentityMap { + template + KOKKOS_INLINE_FUNCTION ValArray value(ValArray const &ref) const { + return ref; + } + + /** + * @brief computes the physical gradient + * + * @details + * With J[i][j] = dx_j/dxi_i the chain rule gives grad_xi = J grad_x, so + * grad_x = J^{-1} grad_xi + */ + template + KOKKOS_INLINE_FUNCTION GradArray gradient(GradArray const &refGrad, + JacobianArray const &J) const { + // sqrt not constexpr until C++26. Is there a more sane way to do this? + constexpr size_t Dim = JacobianArray::size() == 1 ? 1 + : JacobianArray::size() == 4 ? 2 + : 3; + static_assert(Dim * Dim == JacobianArray::size()); + constexpr auto numComp = GradArray::size() / Dim; + const auto Jinv = invert(J); + GradArray g; + for (size_t i = 0; i < Dim; ++i) { + for (size_t ci = 0; ci < numComp; ++ci) { + Real sum = 0; + for (size_t j = 0; j < Dim; ++j) + sum += Jinv[i * Dim + j] * refGrad[j * numComp + ci]; + g[i * numComp + ci] = sum; + } + } + return g; + } +}; + +/** + * @brief Element type used for non-isoparametric, and other elements that require mapping + * + * @details + * The field element supplies the basis and the stored dofs. The geometry + * element supplies the map from the reference element to the physical + * element. I + * + * Evaluation: + * @code + * const auto c = elm.coefficients(ent); // once per element + * for (pt : points) elm.valueFromCoefficients(c, xi); // once per point + * @endcode + * + * @tparam FieldElem a FieldElement over the field being evaluated + * @tparam GeomElem a FieldElement over the coordinate field + * @tparam DofTransformation see IdentityDofTransformation + * @tparam MapType see IdentityMap + */ +template +struct MappedElement { + static_assert(FieldElem::MeshEntDim == GeomElem::MeshEntDim, + "the field and geometry elements must have the same " + "parametric dimension"); + static_assert(FieldElem::getTopology()[0] == GeomElem::getTopology()[0], + "the field and geometry elements must have the same topology"); + + FieldElem field; + GeomElem geom; + DofTransformation dofTransformation; + MapType mapType; + const size_t numMeshEnts; + + static constexpr size_t MeshEntDim = FieldElem::MeshEntDim; + static constexpr size_t NumComponents = FieldElem::NumComponents; + using ValArray = typename FieldElem::ValArray; + + /** @brief the topology of the reference element */ + static constexpr KOKKOS_FUNCTION auto getTopology() { + return GeomElem::getTopology(); + } + + /** @brief physical gradients + * [direction * NumComponents + component] */ + using GradArray = typename FieldElem::LocalGradArray; + /** @brief the Jacobian, J[i][j] = dx_j/dxi_i */ + using JacobianArray = typename GeomElem::JacobianArray; + + /** @brief storage for per-element data*/ + struct Coefficients { + /** coefficients on the reference basis, from the DofTransformation */ + decltype(std::declval().coefficients( + std::declval(), std::declval(), 0)) field; + /** the geometry element's node coordinates */ + typename GeomElem::NodeArray geometry; + }; + + MappedElement( + FieldElem const &fieldIn, GeomElem const &geomIn, + DofTransformation const &dofTransformationIn = DofTransformation(), + MapType const &mapTypeIn = MapType()) + : field(fieldIn), geom(geomIn), dofTransformation(dofTransformationIn), + mapType(mapTypeIn), numMeshEnts(fieldIn.numMeshEnts) { + if (geomIn.numMeshEnts != fieldIn.numMeshEnts) { + fail("the field and geometry elements have the same number of " + "mesh entities (%zu != %zu)\n", + fieldIn.numMeshEnts, geomIn.numMeshEnts); + } + } + + KOKKOS_INLINE_FUNCTION Coefficients coefficients(int ent) const { + return {dofTransformation.coefficients(field, geom, ent), + geom.getNodeValues(ent)}; + } + + /** @brief the physical value of the field at a parametric coordinate */ + KOKKOS_INLINE_FUNCTION ValArray valueFromCoefficients( + Coefficients const &c, Kokkos::Array localCoord) const { + return mapType.value(field.valueFromCoefficients(c.field, localCoord)); + } + + /** @brief the Jacobian of the geometry at a parametric coordinate */ + KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromCoefficients( + Coefficients const &c, Kokkos::Array localCoord) const { + return geom.jacobianFromNodeValues(c.geometry, localCoord); + } + + /** @brief the physical gradient of the field at a parametric coordinate */ + KOKKOS_INLINE_FUNCTION GradArray gradientFromCoefficients( + Coefficients const &c, Kokkos::Array localCoord) const { + return mapType.gradient( + field.localGradientFromCoefficients(c.field, localCoord), + jacobianFromCoefficients(c, localCoord)); + } + + /** + * @brief evaluate the field at a single parametric coordinate + * + * Prefer coefficients and valueFromCoefficients when evaluating an element + * at more than one point. + */ + KOKKOS_INLINE_FUNCTION ValArray + getValue(int ent, Kokkos::Array localCoord) const { + return valueFromCoefficients(coefficients(ent), localCoord); + } + + /** @brief Jacobians, from the geometry element */ + Kokkos::View getJacobians(Kokkos::View localCoords, + Kokkos::View offsets) const { + return geom.getJacobians(localCoords, offsets); + } + + /** @brief Jacobian determinants, from the geometry element */ + template + Kokkos::View getJacobianDeterminants(Matrices const &J) const { + return geom.getJacobianDeterminants(J); + } +}; + +} // namespace MeshField + +#endif diff --git a/test/testOmegahMappedElement.cpp b/test/testOmegahMappedElement.cpp new file mode 100644 index 0000000..0dc0599 --- /dev/null +++ b/test/testOmegahMappedElement.cpp @@ -0,0 +1,231 @@ +#include "MeshField.hpp" +#include "MeshField_Integrate.hpp" +#include "MeshField_MappedElement.hpp" +#include "MeshField_TestFunctions.hpp" +#include "Omega_h_build.hpp" +#include +#include + +// MappedElement tests on triangles and tetrahedra. The meshes are sheared so +// the Jacobians are full matrices. +// +// 1. P2 field on P1 geometry +// 2. finite differences on curved P2 geometry +// 3. measure/Integrator with curved geometry + +using ExecutionSpace = Kokkos::DefaultExecutionSpace; +template +using OmegahMeshField = MeshField::OmegahMeshField; + +constexpr int numPts = 3; + +// points used for testing element values +template KOKKOS_INLINE_FUNCTION Vec point(int pt) { + const MeshField::Real pts[numPts][3] = { + {0.25, 0.25, 0.25}, {0.1, 0.5, 0.15}, {0.4, 0.1, 0.3}}; + Vec xi; + for (size_t d = 0; d < dim; ++d) + xi[d] = pts[pt][d]; + return xi; +} + +// set component c of the field element elm to (c+1)*func at its nodes, with +// the node positions from the coordinate element coordElm +template +void interpolate(FieldElem &elm, CoordElem const &coordElm, Func const &func) { + constexpr size_t numComp = FieldElem::NumComponents; + MeshField::setDofByInterpolation( + elm, KOKKOS_LAMBDA(const int ent, Vec const &xi) { + const auto f = func(coordElm.getValue(ent, xi)); + Kokkos::Array val; + for (size_t c = 0; c < numComp; ++c) + val[c] = (c + 1) * f; + return val; + }); +} + +// 1. P2 field on P1 geometry reproduces a quadratic exactly +template +bool testP2FieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { + QuadraticFunction func; + auto cf1Wc = omf.template CreateLagrangeCoordinateField<1>(Shear()); + constexpr size_t numComp = 2; + auto fieldWc = + omf.template CreateLagrangeField(); + auto geomElm = omf.CreateLagrangeElement(cf1Wc.field); + auto fieldElm = omf.CreateLagrangeElement(fieldWc.field); + interpolate(fieldElm, geomElm, func); + + auto elm = omf.CreateLagrangeElement(fieldWc.field, geomElm); + auto coordElm = omf.CreateLagrangeElement(cf1Wc.field, geomElm); + + MeshField::LO errors = 0; + Kokkos::parallel_reduce( + "checkP2FieldOnP1Geometry", mesh.nelems(), + KOKKOS_LAMBDA(const int ent, MeshField::LO &lerr) { + const auto c = elm.coefficients(ent); + const auto cx = coordElm.coefficients(ent); + for (int pt = 0; pt < numPts; ++pt) { + const auto xi = point(pt); + const auto pos = coordElm.valueFromCoefficients(cx, xi); + Vec x; + for (size_t d = 0; d < dim; ++d) + x[d] = pos[d]; + const auto v = elm.valueFromCoefficients(c, xi); + const auto g = elm.gradientFromCoefficients(c, xi); + const auto gf = func.grad(x); + for (size_t ci = 0; ci < numComp; ++ci) { + if (!(Kokkos::fabs(v[ci] - (ci + 1) * func(x)) <= + MeshField::Epsilon)) + lerr++; + for (size_t d = 0; d < dim; ++d) + if (!(Kokkos::fabs(g[d * numComp + ci] - (ci + 1) * gf[d]) <= + MeshField::Epsilon)) + lerr++; + } + } + }, + errors); + std::cout << "testP2FieldOnP1Geometry(" << dim + << "d): " << (errors ? "FAIL" : "pass") << "\n"; + return errors == 0; +} + +// 2. physical gradients on curved geometry agree with central differences +template +bool testCurvedGradientFiniteDifference(Omega_h::Mesh &mesh, + OmegahMeshField &omf) { + QuadraticFunction func; + // sheared, then curved by a smooth offset + const MeshField::Real s = 0.03; + auto cfWc = omf.template CreateLagrangeCoordinateField<2>( + KOKKOS_LAMBDA(Vec x) { + x = Shear()(x); + auto y = x; + for (size_t d = 0; d < dim; ++d) + y[d] += s * Kokkos::sin(3.0 * x[(d + 1) % dim] + d + 1); + return y; + }); + auto fieldWc = + omf.template CreateLagrangeField(); + auto geomElm = omf.CreateLagrangeElement(cfWc.field); + auto fieldElm = omf.CreateLagrangeElement(fieldWc.field); + interpolate(fieldElm, geomElm, func); + + auto elm = omf.CreateLagrangeElement(fieldWc.field, geomElm); + auto coordElm = omf.CreateLagrangeElement(cfWc.field, geomElm); + + const MeshField::Real h = 1e-5; + MeshField::LO errors = 0; + Kokkos::parallel_reduce( + "checkFiniteDifference", mesh.nelems(), + KOKKOS_LAMBDA(const int ent, MeshField::LO &lerr) { + const auto c = elm.coefficients(ent); + const auto cx = coordElm.coefficients(ent); + for (int pt = 0; pt < numPts; ++pt) { + const auto xi = point(pt); + const auto g = elm.gradientFromCoefficients(c, xi); + for (size_t k = 0; k < dim; ++k) { + auto xp = xi; + auto xm = xi; + xp[k] += h; + xm[k] -= h; + const auto dv = (elm.valueFromCoefficients(c, xp)[0] - + elm.valueFromCoefficients(c, xm)[0]) / + (2 * h); + const auto posP = coordElm.valueFromCoefficients(cx, xp); + const auto posM = coordElm.valueFromCoefficients(cx, xm); + MeshField::Real dvFromGrad = 0; + for (size_t j = 0; j < dim; ++j) + dvFromGrad += (posP[j] - posM[j]) / (2 * h) * g[j]; + if (!(Kokkos::fabs(dvFromGrad - dv) <= 1e-7)) + lerr++; + } + } + }, + errors); + std::cout << "testCurvedGradientFiniteDifference(" << dim << "d, P" + << FieldOrder << " field): " << (errors ? "FAIL" : "pass") << "\n"; + return errors == 0; +} + +// compute the measure (3D volume, 2D area, 1D length) +class MeasureIntegrator : public MeshField::Integrator { +public: + MeasureIntegrator() : MeshField::Integrator(2) {} + void atPoints(Kokkos::View, + Kokkos::View w, + Kokkos::View dV) override { + MeshField::Real sum = 0; + Kokkos::parallel_reduce( + "integrateMeasure", w.extent(0), + KOKKOS_LAMBDA(const int i, MeshField::Real &lsum) { + lsum += w(i) * dV(i); + }, + sum); + measure += sum; + } + MeshField::Real measure = 0; +}; + +// 3. measure of P2 geometry through an Integrator over a MappedElement +template +bool testMeasure(OmegahMeshField &omf, CoordinateMap const &coordMap, + MeshField::Real expected) { + auto cfWc = omf.template CreateLagrangeCoordinateField<2>(coordMap); + auto fieldWc = omf.template CreateLagrangeField(); + auto elm = omf.CreateLagrangeElement(fieldWc.field, + omf.CreateLagrangeElement(cfWc.field)); + MeasureIntegrator integrator; + integrator.process(elm); + const bool pass = + std::abs(integrator.measure - expected) <= MeshField::Epsilon; + std::cout << "testMeasure(" << dim << "d): " << integrator.measure + << " expected " << expected << (pass ? " pass" : " FAIL") << "\n"; + return pass; +} + +template int runTests(Omega_h::Mesh &mesh) { + OmegahMeshField omf(mesh); + int failed = 0; + failed += !testP2FieldOnP1Geometry(mesh, omf); + failed += !testCurvedGradientFiniteDifference(mesh, omf); + failed += !testCurvedGradientFiniteDifference(mesh, omf); + if constexpr (dim == 2) { + // the unit square with its top bent into y = 1 + delta * x(1 - x); the + // top edges are parabolas, so P2 reproduces the boundary and the area + const MeshField::Real delta = 0.1; + failed += !testMeasure( + omf, + KOKKOS_LAMBDA(Vec x) { + x[1] *= 1 + delta * x[0] * (1 - x[0]); + return x; + }, + 1.0 + delta / 6.0); + } else { + failed += !testMeasure(omf, Shear(), 6.0); // 1x2x3 box + } + return failed; +} + +int main(int argc, char **argv) { + Kokkos::initialize(argc, argv); + int failed = 0; + { + auto lib = Omega_h::Library(&argc, &argv); + MeshField::Debug = true; + auto tris = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 1.0, 0.0, + 3, 3, 0); + auto tets = Omega_h::build_box(lib.world(), OMEGA_H_SIMPLEX, 1.0, 2.0, 3.0, + 2, 2, 2); + failed += runTests<2>(tris); + failed += runTests<3>(tets); + } + Kokkos::finalize(); + if (failed) { + std::cerr << failed << " test(s) failed\n"; + return 1; + } + std::cout << "all tests passed\n"; + return 0; +} From 238622edec6ec1b840d19265344b7e911c631fa8 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 26 Sep 2026 22:29:26 -0400 Subject: [PATCH 10/19] support affine and isoparametric MappedElement In the case of affine (straight sided) and isoparametric elements we will want to take compile-time optimized code paths so the performance of a MappedElement is similar to a base FieldElement --- src/MeshField_Element.hpp | 33 ++++++++++++++++++++-- src/MeshField_MappedElement.hpp | 47 ++++++++++++++++++++++++++------ test/testOmegahMappedElement.cpp | 24 +++++++++------- 3 files changed, 84 insertions(+), 20 deletions(-) diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index 341b7b4..916b971 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -274,7 +274,22 @@ struct FieldElement { template KOKKOS_INLINE_FUNCTION LocalGradArray localGradientFromCoefficients( Coefficients const &c, Kokkos::Array localCoord) const { - const auto shapeGrads = shapeFn.getLocalGradients(localCoord); + return localGradientFromShapeGradients( + c, shapeFn.getLocalGradients(localCoord)); + } + + /** + * @brief + * localGradientFromCoefficients with the shape function gradients already + * evaluated, so they can be shared with a geometry that uses the same shape + * + * @param c coefficients from coefficients(ent) + * @param shapeGrads shapeFn.getLocalGradients at the parametric coordinate + * @return parametric gradients, indexed [direction * NumComponents + component] + */ + template + KOKKOS_INLINE_FUNCTION LocalGradArray localGradientFromShapeGradients( + Coefficients const &c, ShapeGradients const &shapeGrads) const { LocalGradArray g; for (size_t i = 0; i < MeshEntDim * NumComponents; ++i) g[i] = 0; @@ -341,7 +356,21 @@ struct FieldElement { */ KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromNodeValues( NodeArray const &x, Kokkos::Array localCoord) const { - const auto grads = shapeFn.getLocalGradients(localCoord); + return jacobianFromShapeGradients(x, shapeFn.getLocalGradients(localCoord)); + } + + /** + * @brief + * jacobianFromNodeValues with the shape function gradients already + * evaluated, so they can be shared with a field that uses the same shape + * + * @param x node coordinates from getNodeValues + * @param grads shapeFn.getLocalGradients at the parametric coordinate + * @return the Jacobian, row-major with J[i][j] = dx_j/dxi_i + */ + template + KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromShapeGradients( + NodeArray const &x, ShapeGradients const &grads) const { JacobianArray J; for (size_t i = 0; i < MeshEntDim * MeshEntDim; ++i) J[i] = 0; diff --git a/src/MeshField_MappedElement.hpp b/src/MeshField_MappedElement.hpp index 8c5477e..2e9ef73 100644 --- a/src/MeshField_MappedElement.hpp +++ b/src/MeshField_MappedElement.hpp @@ -6,6 +6,7 @@ #include #include #include +#include #include // std::declval /** @@ -121,13 +122,25 @@ struct MappedElement { /** @brief the Jacobian, J[i][j] = dx_j/dxi_i */ using JacobianArray = typename GeomElem::JacobianArray; + /** @brief the geometry map is affine, so the Jacobian is constant over + * each element */ + static constexpr bool IsAffine = decltype(GeomElem::shapeFn)::Order == 1; + /** @brief the field and geometry use the same shape functions and dof map */ + static constexpr bool IsIsoparametric = + std::is_same_v && + std::is_same_v; + /** @brief storage for per-element data*/ struct Coefficients { /** coefficients on the reference basis, from the DofTransformation */ decltype(std::declval().coefficients( std::declval(), std::declval(), 0)) field; - /** the geometry element's node coordinates */ - typename GeomElem::NodeArray geometry; + /** the Jacobian for affine geometry, which needs fewer values than the + * node coordinates; otherwise the node coordinates */ + std::conditional_t + geometry; }; MappedElement( @@ -144,8 +157,14 @@ struct MappedElement { } KOKKOS_INLINE_FUNCTION Coefficients coefficients(int ent) const { - return {dofTransformation.coefficients(field, geom, ent), - geom.getNodeValues(ent)}; + if constexpr (IsAffine) { + return {dofTransformation.coefficients(field, geom, ent), + geom.jacobianFromNodeValues(geom.getNodeValues(ent), + Kokkos::Array{})}; + } else { + return {dofTransformation.coefficients(field, geom, ent), + geom.getNodeValues(ent)}; + } } /** @brief the physical value of the field at a parametric coordinate */ @@ -157,15 +176,27 @@ struct MappedElement { /** @brief the Jacobian of the geometry at a parametric coordinate */ KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromCoefficients( Coefficients const &c, Kokkos::Array localCoord) const { - return geom.jacobianFromNodeValues(c.geometry, localCoord); + if constexpr (IsAffine) { + return c.geometry; + } else { + return geom.jacobianFromNodeValues(c.geometry, localCoord); + } } /** @brief the physical gradient of the field at a parametric coordinate */ KOKKOS_INLINE_FUNCTION GradArray gradientFromCoefficients( Coefficients const &c, Kokkos::Array localCoord) const { - return mapType.gradient( - field.localGradientFromCoefficients(c.field, localCoord), - jacobianFromCoefficients(c, localCoord)); + if constexpr (IsIsoparametric && !IsAffine) { + // the field and the geometry share the shape function gradients + const auto grads = field.shapeFn.getLocalGradients(localCoord); + return mapType.gradient( + field.localGradientFromShapeGradients(c.field, grads), + geom.jacobianFromShapeGradients(c.geometry, grads)); + } else { + return mapType.gradient( + field.localGradientFromCoefficients(c.field, localCoord), + jacobianFromCoefficients(c, localCoord)); + } } /** diff --git a/test/testOmegahMappedElement.cpp b/test/testOmegahMappedElement.cpp index 0dc0599..0132c41 100644 --- a/test/testOmegahMappedElement.cpp +++ b/test/testOmegahMappedElement.cpp @@ -5,11 +5,12 @@ #include "Omega_h_build.hpp" #include #include +#include // MappedElement tests on triangles and tetrahedra. The meshes are sheared so // the Jacobians are full matrices. // -// 1. P2 field on P1 geometry +// 1. P1 and P2 fields on P1 geometry // 2. finite differences on curved P2 geometry // 3. measure/Integrator with curved geometry @@ -44,14 +45,16 @@ void interpolate(FieldElem &elm, CoordElem const &coordElm, Func const &func) { }); } -// 1. P2 field on P1 geometry reproduces a quadratic exactly -template -bool testP2FieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { - QuadraticFunction func; +// 1. a field on P1 geometry reproduces a polynomial of its order exactly +template +bool testFieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { + std::conditional_t, + QuadraticFunction> + func; auto cf1Wc = omf.template CreateLagrangeCoordinateField<1>(Shear()); constexpr size_t numComp = 2; auto fieldWc = - omf.template CreateLagrangeField(); + omf.template CreateLagrangeField(); auto geomElm = omf.CreateLagrangeElement(cf1Wc.field); auto fieldElm = omf.CreateLagrangeElement(fieldWc.field); interpolate(fieldElm, geomElm, func); @@ -61,7 +64,7 @@ bool testP2FieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { MeshField::LO errors = 0; Kokkos::parallel_reduce( - "checkP2FieldOnP1Geometry", mesh.nelems(), + "checkFieldOnP1Geometry", mesh.nelems(), KOKKOS_LAMBDA(const int ent, MeshField::LO &lerr) { const auto c = elm.coefficients(ent); const auto cx = coordElm.coefficients(ent); @@ -86,8 +89,8 @@ bool testP2FieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { } }, errors); - std::cout << "testP2FieldOnP1Geometry(" << dim - << "d): " << (errors ? "FAIL" : "pass") << "\n"; + std::cout << "testFieldOnP1Geometry(" << dim << "d, P" << FieldOrder + << " field): " << (errors ? "FAIL" : "pass") << "\n"; return errors == 0; } @@ -188,7 +191,8 @@ bool testMeasure(OmegahMeshField &omf, CoordinateMap const &coordMap, template int runTests(Omega_h::Mesh &mesh) { OmegahMeshField omf(mesh); int failed = 0; - failed += !testP2FieldOnP1Geometry(mesh, omf); + failed += !testFieldOnP1Geometry(mesh, omf); + failed += !testFieldOnP1Geometry(mesh, omf); failed += !testCurvedGradientFiniteDifference(mesh, omf); failed += !testCurvedGradientFiniteDifference(mesh, omf); if constexpr (dim == 2) { From 7b9ea57cbaca71dee4f6145a8cd747d67054c3a2 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sat, 26 Sep 2026 23:02:24 -0400 Subject: [PATCH 11/19] gather MappedElement geometry --- src/MeshField_MappedElement.hpp | 62 +++++++++++++++++--------------- test/testOmegahMappedElement.cpp | 6 ++-- 2 files changed, 37 insertions(+), 31 deletions(-) diff --git a/src/MeshField_MappedElement.hpp b/src/MeshField_MappedElement.hpp index 2e9ef73..9cff8e5 100644 --- a/src/MeshField_MappedElement.hpp +++ b/src/MeshField_MappedElement.hpp @@ -82,8 +82,12 @@ struct IdentityMap { * * Evaluation: * @code - * const auto c = elm.coefficients(ent); // once per element - * for (pt : points) elm.valueFromCoefficients(c, xi); // once per point + * const auto c = elm.coefficients(ent); // once per element + * const auto g = elm.geometry(ent); // if needed + * for (pt : points) { + * elm.valueFromCoefficients(c, xi); // values + * elm.gradientFromCoefficients(c, g, xi); // gradients + * } * @endcode * * @tparam FieldElem a FieldElement over the field being evaluated @@ -132,16 +136,11 @@ struct MappedElement { std::is_same_v; - /** @brief storage for per-element data*/ - struct Coefficients { - /** coefficients on the reference basis, from the DofTransformation */ - decltype(std::declval().coefficients( - std::declval(), std::declval(), 0)) field; - /** the Jacobian for affine geometry, which needs fewer values than the - * node coordinates; otherwise the node coordinates */ - std::conditional_t - geometry; - }; + /** @brief the field's coefficients on the reference basis */ + using Coefficients = decltype(std::declval().coefficients( + std::declval(), std::declval(), 0)); + using Geometry = + std::conditional_t; MappedElement( FieldElem const &fieldIn, GeomElem const &geomIn, @@ -157,45 +156,50 @@ struct MappedElement { } KOKKOS_INLINE_FUNCTION Coefficients coefficients(int ent) const { + return dofTransformation.coefficients(field, geom, ent); + } + + /** + * @brief the geometry of one element + */ + KOKKOS_INLINE_FUNCTION Geometry geometry(int ent) const { if constexpr (IsAffine) { - return {dofTransformation.coefficients(field, geom, ent), - geom.jacobianFromNodeValues(geom.getNodeValues(ent), - Kokkos::Array{})}; + return geom.jacobianFromNodeValues(geom.getNodeValues(ent), + Kokkos::Array{}); } else { - return {dofTransformation.coefficients(field, geom, ent), - geom.getNodeValues(ent)}; + return geom.getNodeValues(ent); } } /** @brief the physical value of the field at a parametric coordinate */ KOKKOS_INLINE_FUNCTION ValArray valueFromCoefficients( Coefficients const &c, Kokkos::Array localCoord) const { - return mapType.value(field.valueFromCoefficients(c.field, localCoord)); + return mapType.value(field.valueFromCoefficients(c, localCoord)); } /** @brief the Jacobian of the geometry at a parametric coordinate */ - KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromCoefficients( - Coefficients const &c, Kokkos::Array localCoord) const { + KOKKOS_INLINE_FUNCTION JacobianArray jacobianFromGeometry( + Geometry const &g, Kokkos::Array localCoord) const { if constexpr (IsAffine) { - return c.geometry; + return g; } else { - return geom.jacobianFromNodeValues(c.geometry, localCoord); + return geom.jacobianFromNodeValues(g, localCoord); } } /** @brief the physical gradient of the field at a parametric coordinate */ - KOKKOS_INLINE_FUNCTION GradArray gradientFromCoefficients( - Coefficients const &c, Kokkos::Array localCoord) const { + KOKKOS_INLINE_FUNCTION GradArray + gradientFromCoefficients(Coefficients const &c, Geometry const &g, + Kokkos::Array localCoord) const { if constexpr (IsIsoparametric && !IsAffine) { // the field and the geometry share the shape function gradients const auto grads = field.shapeFn.getLocalGradients(localCoord); - return mapType.gradient( - field.localGradientFromShapeGradients(c.field, grads), - geom.jacobianFromShapeGradients(c.geometry, grads)); + return mapType.gradient(field.localGradientFromShapeGradients(c, grads), + geom.jacobianFromShapeGradients(g, grads)); } else { return mapType.gradient( - field.localGradientFromCoefficients(c.field, localCoord), - jacobianFromCoefficients(c, localCoord)); + field.localGradientFromCoefficients(c, localCoord), + jacobianFromGeometry(g, localCoord)); } } diff --git a/test/testOmegahMappedElement.cpp b/test/testOmegahMappedElement.cpp index 0132c41..3761e86 100644 --- a/test/testOmegahMappedElement.cpp +++ b/test/testOmegahMappedElement.cpp @@ -67,6 +67,7 @@ bool testFieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { "checkFieldOnP1Geometry", mesh.nelems(), KOKKOS_LAMBDA(const int ent, MeshField::LO &lerr) { const auto c = elm.coefficients(ent); + const auto geo = elm.geometry(ent); const auto cx = coordElm.coefficients(ent); for (int pt = 0; pt < numPts; ++pt) { const auto xi = point(pt); @@ -75,7 +76,7 @@ bool testFieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { for (size_t d = 0; d < dim; ++d) x[d] = pos[d]; const auto v = elm.valueFromCoefficients(c, xi); - const auto g = elm.gradientFromCoefficients(c, xi); + const auto g = elm.gradientFromCoefficients(c, geo, xi); const auto gf = func.grad(x); for (size_t ci = 0; ci < numComp; ++ci) { if (!(Kokkos::fabs(v[ci] - (ci + 1) * func(x)) <= @@ -124,10 +125,11 @@ bool testCurvedGradientFiniteDifference(Omega_h::Mesh &mesh, "checkFiniteDifference", mesh.nelems(), KOKKOS_LAMBDA(const int ent, MeshField::LO &lerr) { const auto c = elm.coefficients(ent); + const auto geo = elm.geometry(ent); const auto cx = coordElm.coefficients(ent); for (int pt = 0; pt < numPts; ++pt) { const auto xi = point(pt); - const auto g = elm.gradientFromCoefficients(c, xi); + const auto g = elm.gradientFromCoefficients(c, geo, xi); for (size_t k = 0; k < dim; ++k) { auto xp = xi; auto xm = xi; From e6d561854922581ffc6351f9dfd86a271bad3893 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sun, 11 Oct 2026 14:45:41 -0400 Subject: [PATCH 12/19] add MappedElement overloads of setDofByInterpolation --- src/MeshField.hpp | 22 ++-------------- src/MeshField_Element.hpp | 3 +++ src/MeshField_MappedElement.hpp | 37 +++++++++++++++++++++++++++ test/MeshField_TestFunctions.hpp | 28 +++++++++++++++++++++ test/testInterpolate.cpp | 13 ++++------ test/testOmegahMappedElement.cpp | 43 ++++++++++---------------------- 6 files changed, 88 insertions(+), 58 deletions(-) diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 5202ff2..e67421b 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -54,21 +54,6 @@ struct SetCoordFieldFunctor const CoordinateMap map_; }; -template -struct InterpolateElementFunctor -{ - InterpolateElementFunctor(const Element& element, const CoordinateMap& map) - : element_(element), map_(map) {} - - KOKKOS_INLINE_FUNCTION - auto operator()(const int ent, Kokkos::Array const& xi) const { - return map_(element_.getValue(ent, xi)); - } - - const Element element_; - const CoordinateMap map_; -}; - } // anonymous namespace namespace MeshField { @@ -365,11 +350,8 @@ class OmegahMeshField { "setCoordField"); } else { const auto linear = CreateLagrangeElement(coordField.field); - auto target = CreateLagrangeElement(field); - MeshField::setDofByInterpolation( - target, - InterpolateElementFunctor( - linear, map)); + auto target = CreateLagrangeElement(field, linear); + MeshField::setDofByInterpolation(target, map); } return fieldWithCtrlr; } diff --git a/src/MeshField_Element.hpp b/src/MeshField_Element.hpp index 916b971..6adc30c 100644 --- a/src/MeshField_Element.hpp +++ b/src/MeshField_Element.hpp @@ -809,6 +809,9 @@ DofHolderOwningEntities getDofHolderOwningEntities(FieldElement const &fes) { * @brief * evaluates a source function at the target field element dof holder locations * after this function runs, the target's field values will modified + * + * @details + * Reference space interpolation. To interpolate in physical space see MappedElement overload. */ template void setDofByInterpolation(FieldElement &target, Source const &source) { diff --git a/src/MeshField_MappedElement.hpp b/src/MeshField_MappedElement.hpp index 9cff8e5..93d1d9f 100644 --- a/src/MeshField_MappedElement.hpp +++ b/src/MeshField_MappedElement.hpp @@ -227,6 +227,43 @@ struct MappedElement { } }; +/** + * @brief + * evaluates a source function at the dof holder locations of the target's + * field, in reference space + */ +template + requires std::is_invocable_v const &> +void setDofByInterpolation( + MappedElement &target, + Source const &source) { + setDofByInterpolation(target.field, source); +} + +/** + * @brief + * evaluates a function at the dof holder locations of the target's field, in + * physical space + * + */ +template + requires std::is_invocable_v const &> +void setDofByInterpolation( + MappedElement &target, + Function const &f) { + const auto geom = target.geom; + setDofByInterpolation( + target.field, + KOKKOS_LAMBDA(const int ent, + Kokkos::Array const &xi) { + return f(geom.getValue(ent, xi)); + }); +} + } // namespace MeshField #endif diff --git a/test/MeshField_TestFunctions.hpp b/test/MeshField_TestFunctions.hpp index e0997a3..5e6d8f2 100644 --- a/test/MeshField_TestFunctions.hpp +++ b/test/MeshField_TestFunctions.hpp @@ -50,6 +50,34 @@ template struct QuadraticFunction { } }; +/** + * @brief the value of a numComp component function whose component c is + * (c+1)*f, which checks the component ordering + */ +template +KOKKOS_INLINE_FUNCTION Kokkos::Array +valueScaledByComponent(Function const &f, Vec const &x) { + Kokkos::Array val; + for (size_t c = 0; c < numComp; ++c) + val[c] = (c + 1) * f(x); + return val; +} + +/** + * @brief the gradient of the function in valueScaledByComponent, indexed + * [direction * numComp + component]; f must provide grad + */ +template +KOKKOS_INLINE_FUNCTION Kokkos::Array +gradScaledByComponent(Function const &f, Vec const &x) { + const auto gf = f.grad(x); + Kokkos::Array g; + for (size_t d = 0; d < dim; ++d) + for (size_t c = 0; c < numComp; ++c) + g[d * numComp + c] = (c + 1) * gf[d]; + return g; +} + /** * @brief CoordinateMap for an affine shear that makes every Jacobian a full * matrix: x += 0.3y + 0.1z, y += 0.2z diff --git a/test/testInterpolate.cpp b/test/testInterpolate.cpp index 71233f1..afb35df 100644 --- a/test/testInterpolate.cpp +++ b/test/testInterpolate.cpp @@ -53,13 +53,10 @@ bool testAnalytic( const auto f = KOKKOS_LAMBDA(Kokkos::Array const &x) { return order == 1 ? LinearFunction()(x) : QuadraticFunction()(x); }; + MeshField::MappedElement mappedElm(elm, coordElm); MeshField::setDofByInterpolation( - elm, KOKKOS_LAMBDA(const int ent, Kokkos::Array const &xi) { - const auto x = coordElm.getValue(ent, xi); - Kokkos::Array val; - for (size_t c = 0; c < numComp; ++c) - val[c] = (c + 1) * f(x); - return val; + mappedElm, KOKKOS_LAMBDA(Kokkos::Array const &x) { + return valueScaledByComponent(f, x); }); const auto lc = offNodeCoords(mesh.nelems()); @@ -72,9 +69,9 @@ bool testAnalytic( Kokkos::Array xp; for (size_t d = 0; d < dim; ++d) xp[d] = x(pt, d); + const auto expected = valueScaledByComponent(f, xp); for (size_t c = 0; c < numComp; ++c) - if (Kokkos::fabs(values(pt, c) - (c + 1) * f(xp)) > - MeshField::Epsilon) + if (Kokkos::fabs(values(pt, c) - expected[c]) > MeshField::Epsilon) ++lerrors; }, numErrors); diff --git a/test/testOmegahMappedElement.cpp b/test/testOmegahMappedElement.cpp index 3761e86..54871ff 100644 --- a/test/testOmegahMappedElement.cpp +++ b/test/testOmegahMappedElement.cpp @@ -30,21 +30,6 @@ template KOKKOS_INLINE_FUNCTION Vec point(int pt) { return xi; } -// set component c of the field element elm to (c+1)*func at its nodes, with -// the node positions from the coordinate element coordElm -template -void interpolate(FieldElem &elm, CoordElem const &coordElm, Func const &func) { - constexpr size_t numComp = FieldElem::NumComponents; - MeshField::setDofByInterpolation( - elm, KOKKOS_LAMBDA(const int ent, Vec const &xi) { - const auto f = func(coordElm.getValue(ent, xi)); - Kokkos::Array val; - for (size_t c = 0; c < numComp; ++c) - val[c] = (c + 1) * f; - return val; - }); -} - // 1. a field on P1 geometry reproduces a polynomial of its order exactly template bool testFieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { @@ -56,10 +41,10 @@ bool testFieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { auto fieldWc = omf.template CreateLagrangeField(); auto geomElm = omf.CreateLagrangeElement(cf1Wc.field); - auto fieldElm = omf.CreateLagrangeElement(fieldWc.field); - interpolate(fieldElm, geomElm, func); - auto elm = omf.CreateLagrangeElement(fieldWc.field, geomElm); + MeshField::setDofByInterpolation(elm, KOKKOS_LAMBDA(Vec const &x) { + return valueScaledByComponent(func, x); + }); auto coordElm = omf.CreateLagrangeElement(cf1Wc.field, geomElm); MeshField::LO errors = 0; @@ -77,16 +62,14 @@ bool testFieldOnP1Geometry(Omega_h::Mesh &mesh, OmegahMeshField &omf) { x[d] = pos[d]; const auto v = elm.valueFromCoefficients(c, xi); const auto g = elm.gradientFromCoefficients(c, geo, xi); - const auto gf = func.grad(x); - for (size_t ci = 0; ci < numComp; ++ci) { - if (!(Kokkos::fabs(v[ci] - (ci + 1) * func(x)) <= - MeshField::Epsilon)) + const auto expected = valueScaledByComponent(func, x); + const auto expectedGrad = gradScaledByComponent(func, x); + for (size_t i = 0; i < numComp; ++i) + if (!(Kokkos::fabs(v[i] - expected[i]) <= MeshField::Epsilon)) + lerr++; + for (size_t i = 0; i < dim * numComp; ++i) + if (!(Kokkos::fabs(g[i] - expectedGrad[i]) <= MeshField::Epsilon)) lerr++; - for (size_t d = 0; d < dim; ++d) - if (!(Kokkos::fabs(g[d * numComp + ci] - (ci + 1) * gf[d]) <= - MeshField::Epsilon)) - lerr++; - } } }, errors); @@ -113,10 +96,10 @@ bool testCurvedGradientFiniteDifference(Omega_h::Mesh &mesh, auto fieldWc = omf.template CreateLagrangeField(); auto geomElm = omf.CreateLagrangeElement(cfWc.field); - auto fieldElm = omf.CreateLagrangeElement(fieldWc.field); - interpolate(fieldElm, geomElm, func); - auto elm = omf.CreateLagrangeElement(fieldWc.field, geomElm); + MeshField::setDofByInterpolation(elm, KOKKOS_LAMBDA(Vec const &x) { + return valueScaledByComponent<1>(func, x); + }); auto coordElm = omf.CreateLagrangeElement(cfWc.field, geomElm); const MeshField::Real h = 1e-5; From a9469ca29317f999b13b455014e7d2e5f2563979 Mon Sep 17 00:00:00 2001 From: Jacob Merson Date: Sun, 11 Oct 2026 16:33:36 -0400 Subject: [PATCH 13/19] Add TestHarness to test with all controllers/mesh/elements The TestHarness runs tests with all available controllers and mesh types. This will improve generalizability as we add additional mesh types e.g., mixed and new element types e.g., quads --- test/MeshField_TestFunctions.hpp | 16 +++ test/MeshField_TestHarness.hpp | 162 +++++++++++++++++++++++++++++++ test/testIntegratorMeasure.cpp | 119 ++++++++++------------- test/testInterpolate.cpp | 98 ++++++------------- test/testOmegahMappedElement.cpp | 130 ++++++------------------- 5 files changed, 292 insertions(+), 233 deletions(-) create mode 100644 test/MeshField_TestHarness.hpp diff --git a/test/MeshField_TestFunctions.hpp b/test/MeshField_TestFunctions.hpp index 5e6d8f2..4559900 100644 --- a/test/MeshField_TestFunctions.hpp +++ b/test/MeshField_TestFunctions.hpp @@ -108,4 +108,20 @@ template struct Bend { } }; +/** + * @brief CoordinateMap that bulges the top of the unit square into + * y = 1 + delta x(1 - x): y *= 1 + delta x(1 - x) + * + * The top edges are parabolas in their edge parameter, so quadratic + * coordinate fields reproduce the boundary, and the area, 1 + delta/6, + * exactly. + */ +template struct Bulge { + MeshField::Real delta = 0.1; + KOKKOS_INLINE_FUNCTION Vec operator()(Vec x) const { + x[1] *= 1 + delta * x[0] * (1 - x[0]); + return x; + } +}; + #endif diff --git a/test/MeshField_TestHarness.hpp b/test/MeshField_TestHarness.hpp new file mode 100644 index 0000000..5a99934 --- /dev/null +++ b/test/MeshField_TestHarness.hpp @@ -0,0 +1,162 @@ +#ifndef MESHFIELD_TESTHARNESS_HPP +#define MESHFIELD_TESTHARNESS_HPP + +/** + * @file MeshField_TestHarness.hpp + * @brief a test runner, reporting, and integrators shared by the tests + * + * A test's main creates a TestRunner and loops over controllers, meshes, and + * types (e.g., element factories) with its forEach functions. The loop + * bodies are template lambdas: + * + * @code + * TestRunner t(argc, argv); + * t.forEachController([&]