From c0bde8a870245269642f3783ee74d66bbbd6ea96 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 6 Aug 2026 14:01:38 -0400 Subject: [PATCH 01/22] docs: adding shape functions [WIP] --- Doxyfile.in | 1 + docs/MeshFields.dox | 4 + docs/addingShapeFunctions.md | 18 -- docs/addingShapeFunctionsForOmegah.md | 377 ++++++++++++++++++++++++++ 4 files changed, 382 insertions(+), 18 deletions(-) delete mode 100644 docs/addingShapeFunctions.md create mode 100644 docs/addingShapeFunctionsForOmegah.md diff --git a/Doxyfile.in b/Doxyfile.in index 245a090d..d6f1bce2 100644 --- a/Doxyfile.in +++ b/Doxyfile.in @@ -919,6 +919,7 @@ WARN_LOGFILE = INPUT = \ @CMAKE_CURRENT_SOURCE_DIR@/docs/MeshFields.dox \ +@CMAKE_CURRENT_SOURCE_DIR@/docs/addingShapeFunctionsForOmegah.md \ @CMAKE_CURRENT_SOURCE_DIR@/src/MeshField_Field.hpp \ @CMAKE_CURRENT_SOURCE_DIR@/src/MeshField_Shape.hpp \ @CMAKE_CURRENT_SOURCE_DIR@/src/MeshField_ShapeField.hpp \ diff --git a/docs/MeshFields.dox b/docs/MeshFields.dox index 0ea3ce2b..cf0566ca 100644 --- a/docs/MeshFields.dox +++ b/docs/MeshFields.dox @@ -19,6 +19,10 @@ * exist per mesh entity. For example, a mesh edge could have multiple * nodes for a high order shape function. * + * ## Developer Guides + * + * - @subpage adding-shape-functions-omegah + * * ## References and Related Work * * - PUMI users guide: https://www.scorec.rpi.edu/pumi/PUMI.pdf diff --git a/docs/addingShapeFunctions.md b/docs/addingShapeFunctions.md deleted file mode 100644 index 8f235012..00000000 --- a/docs/addingShapeFunctions.md +++ /dev/null @@ -1,18 +0,0 @@ -In order to add shape function, you must create a struct in -https://github.com/SCOREC/meshFields/blob/main/src/MeshField_Shape.hpp. - -You must define numNodes, meshEntDim, and Order as const variables as well as the functions -getValues and getLocalGradients. You then need to define a mapping from how you store the -nodes to the ordering of the shape function defined earlier. The Omegah mappings are in -https://github.com/SCOREC/meshFields/blob/main/src/MeshField.hpp#L63. In the Jacobian -tests we also an identity mapping defined as another example: -https://github.com/SCOREC/meshFields/blob/main/test/testElementJacobian2d.cpp#L12. - -In order to use integration with the added shape function it must be added here: -https://github.com/SCOREC/meshFields/blob/main/src/MeshField_Integrate.hpp#L44. - -Finally, when trying to compute an integral using the shape function, a class derived from the -class Integrator (defined here: -https://github.com/SCOREC/meshFields/blob/main/src/MeshField_Integrate.hpp#L204) must be -declared and should implement the atPoints function. An example of this can be found here: -https://github.com/SCOREC/meshFields/blob/main/test/testCountIntegrator.cpp#L41. diff --git a/docs/addingShapeFunctionsForOmegah.md b/docs/addingShapeFunctionsForOmegah.md new file mode 100644 index 00000000..10ddcd04 --- /dev/null +++ b/docs/addingShapeFunctionsForOmegah.md @@ -0,0 +1,377 @@ +# Adding Shape Functions for Omega_h Element Topologies {#adding-shape-functions-omegah} + +This guide walks through every step required to add a new shape function +for an Omega\_h element topology (triangle, tetrahedron, quad, etc.). +Four source locations must be modified in order: + +1. Define the **shape function struct** in `src/MeshField_Shape.hpp` +2. Define the **Omega\_h node mapping struct** in `src/MeshField.hpp` +3. Register the shape + mapping in the **element factory** in `src/MeshField.hpp` +4. Add **integration support** in `src/MeshField_Integrate.hpp` + +An optional fifth step covers writing an `Integrator`-derived class to +perform numerical integration with the new shape. + +--- + +## Background + +### Parametric Coordinates and Node Ordering + +meshFields defines shape functions in a canonical parametric coordinate +system that is independent of Omega\_h. Omega\_h uses its own node +ordering for element vertices and edges, which may differ from the +meshFields canonical ordering. The **node mapping struct** (Step 2) +bridges the two conventions. + +For simplex elements the parametric coordinates are reduced barycentric +coordinates. The redundant coordinate \f$L_0 = 1 - \sum \xi_i\f$ is +omitted: + +| Topology | Parametric coords | Range | +|---------------|-------------------|-----------------| +| Edge (1D) | \f$\xi\f$ | \f$[-1, 1]\f$ | +| Triangle (2D) | \f$(\xi_0,\xi_1)\f$ | \f$[0,1]^2\f$, \f$\xi_0+\xi_1 \le 1\f$ | +| Tet (3D) | \f$(\xi_0,\xi_1,\xi_2)\f$ | \f$[0,1]^3\f$, \f$\sum \xi_i \le 1\f$ | + +### Terminology + +- **Node** – a parametric location on an element that carries a DOF. +- **DOF holder** – the mesh entity (Vertex, Edge, …) that owns a DOF. + For a linear triangle the three nodes coincide with the three vertices. + For a quadratic triangle there are also three mid-edge nodes. +- **Node mapping** – a callable that, given a node index local to an + element, returns the global mesh entity index and topology type that + holds the corresponding DOF. + +--- + +## Step 1 – Define the Shape Function Struct + +Add a new struct to `src/MeshField_Shape.hpp` inside `namespace MeshField`. + +### Required Members + +| Member | Type | Description | +|--------|------|-------------| +| `numNodes` | `static const size_t` | Total nodes per element | +| `meshEntDim` | `static const size_t` | Parametric space dimension | +| `Order` | `constexpr static size_t` | Polynomial order | +| `DofHolders` | `constexpr static Mesh_Topology[]` | Entity types that hold DOFs | +| `getNodeParametricCoords()` | `KOKKOS_INLINE_FUNCTION` | Flat array of node coordinates (length `numNodes * meshEntDim`) | +| `getValues(xi)` | `KOKKOS_INLINE_FUNCTION` | Shape function values at `xi` (length `numNodes`) | +| `getLocalGradients(xi)` | `KOKKOS_INLINE_FUNCTION` | Gradients at `xi` (length `meshEntDim * numNodes`, row-major: \f$[\partial N_0/\partial\xi_0, \partial N_0/\partial\xi_1, \ldots]\f$) | + +Optional (needed by quadratic and higher-order shapes): + +| Member | Type | Description | +|--------|------|-------------| +| `NumDofHolders` | `constexpr static size_t[]` | Count of entities of each type in `DofHolders` | +| `DofsPerHolder` | `constexpr static size_t[]` | DOFs per entity for each type in `DofHolders` | + +### Validating Parametric Coordinates + +Use the helper functions declared at the top of `MeshField_Shape.hpp` to +assert that incoming parametric coordinates are in range. Pass +`MeshField::ParametricCoordTol` as the tolerance. + +```cpp +assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol)); +assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol)); +// for barycentric remainder: +const Real L0 = 1 - xi[0] - xi[1]; +assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); +``` + +### Example Skeleton (Linear Quadrilateral) + +```cpp +struct LinearQuadShape { + static const size_t numNodes = 4; + static const size_t meshEntDim = 2; + constexpr static Mesh_Topology DofHolders[1] = {Vertex}; + constexpr static size_t Order = 1; + + KOKKOS_INLINE_FUNCTION + Kokkos::Array getNodeParametricCoords() const { + // clang-format off + return {-1,-1, // node 0 + 1,-1, // node 1 + 1, 1, // node 2 + -1, 1}; // node 3 + // clang-format on + } + + KOKKOS_INLINE_FUNCTION + Kokkos::Array getValues(Vector2 const &xi) const { + // standard bilinear shape functions + const Real xim = 1.0 - xi[0], xip = 1.0 + xi[0]; + const Real etam = 1.0 - xi[1], etap = 1.0 + xi[1]; + // clang-format off + return {0.25*xim*etam, + 0.25*xip*etam, + 0.25*xip*etap, + 0.25*xim*etap}; + // clang-format on + } + + KOKKOS_INLINE_FUNCTION + Kokkos::Array getLocalGradients(Vector2 const &xi) const { + const Real xim = 1.0 - xi[0], xip = 1.0 + xi[0]; + const Real etam = 1.0 - xi[1], etap = 1.0 + xi[1]; + // clang-format off + return {-0.25*etam, -0.25*xim, // dN0/dxi0, dN0/dxi1 + 0.25*etam, -0.25*xip, // dN1/dxi0, dN1/dxi1 + 0.25*etap, 0.25*xip, // dN2/dxi0, dN2/dxi1 + -0.25*etap, 0.25*xim}; // dN3/dxi0, dN3/dxi1 + // clang-format on + } +}; +``` + +--- + +## Step 2 – Define the Omega\_h Node Mapping Struct + +The mapping struct lives in `src/MeshField.hpp` inside +`namespace MeshField::Omegah`. It translates element-local node indices +(as numbered by the shape function) into the global mesh entity indices +(as numbered by Omega\_h) that hold the corresponding DOFs. + +### Interface Requirements + +| Member | Description | +|--------|-------------| +| Constructor `(Omega_h::Mesh &)` | Cache connectivity arrays from Omega\_h; validate mesh family/dimension | +| `static constexpr getTopology()` | Return a `Kokkos::Array` listing element topologies this mapping applies to | +| `operator()(LO nodeIdx, LO compIdx, LO elem, Mesh_Topology topo)` | Return `ElementToDofHolderMap{node, comp, ent, topo}` | + +`ElementToDofHolderMap` packs four values: `{nodeLocalIdx, componentIdx, globalEntityIdx, entityTopology}`. +For typical shape functions `nodeLocalIdx` is always `0`. + +### Omega\_h Connectivity APIs + +| Query | API call | +|-------|----------| +| Element→vertex connectivity | `mesh.ask_elem_verts()` → flat `LOs` array, stride = `simplex_degree(elemDim, 0)` | +| Element→edge connectivity | `mesh.ask_down(elemDim, 1).ab2b` → flat `LOs` array, stride = `simplex_degree(elemDim, 1)` | +| Element→face connectivity | `mesh.ask_down(elemDim, 2).ab2b` | + +### Vertex-Ordering Correction + +Omega\_h numbers vertices and edges within a simplex differently from the +meshFields canonical ordering used by the shape functions. +The existing linear-triangle and linear-tetrahedron mappings correct for +this with a cyclic rotation: + +```cpp +// For triangles (triDim=2, vtxDim=0): +const auto localVtxIdx = + (Omega_h::simplex_down_template(triDim, vtxDim, nodeIdx, /*ignored=*/-1) + 2) % 3; + +// For tetrahedra (tetDim=3, vtxDim=0): +const auto localVtxIdx = + (Omega_h::simplex_down_template(tetDim, vtxDim, nodeIdx, /*ignored=*/-1) + 3) % 4; +``` + +You must determine the correct rotation offset for your topology by +comparing the Omega\_h canonical vertex/edge ordering (see +`Omega_h_simplex.hpp` and the `simplex_down_template` function) with the +node ordering established by `getNodeParametricCoords()` in your shape +struct. Validate with a unit test that evaluates the shape functions at +the parametric coordinates of each node and checks that the value for +that node is 1.0 and all others are 0.0. + +### Example Skeleton (Linear Quad, vertices only) + +```cpp +struct LinearQuadToVertexField { + Omega_h::LOs quadVerts; + + LinearQuadToVertexField(Omega_h::Mesh &mesh) + : quadVerts(mesh.ask_elem_verts()) { + if (mesh.dim() != 2 || mesh.family() != OMEGA_H_HYPERCUBE) + MeshField::fail("Mesh must be 2D hypercube (quads)\n"); + } + + static constexpr KOKKOS_FUNCTION Kokkos::Array + getTopology() { + return {MeshField::Quad}; + } + + KOKKOS_FUNCTION MeshField::ElementToDofHolderMap + operator()(MeshField::LO nodeIdx, MeshField::LO compIdx, + MeshField::LO elem, MeshField::Mesh_Topology topo) const { + assert(topo == MeshField::Quad); + const auto quadDim = 2, vtxDim = 0; + // Determine the rotation offset by comparing Omega_h vertex ordering + // with your shape function's getNodeParametricCoords() ordering. + const auto localVtxIdx = + (Omega_h::simplex_down_template(quadDim, vtxDim, nodeIdx, -1) + OFFSET) % 4; + const auto stride = Omega_h::simplex_degree(quadDim, vtxDim); + const MeshField::LO vtx = quadVerts[elem * stride + localVtxIdx]; + return {0, compIdx, vtx, MeshField::Vertex}; + } +}; +``` + +Replace `OFFSET` with the value (0–3) that maps Omega\_h's vertex order +to the meshFields canonical order for your element. See the existing +`LinearTriangleToVertexField` (offset 2) and +`LinearTetrahedronToVertexField` (offset 3) in `src/MeshField.hpp` for +reference. + +--- + +## Step 3 – Register in the Element Factory + +Add a new factory function (or extend an existing one) in +`src/MeshField.hpp` inside `namespace MeshField::Omegah`. +The existing helpers follow the pattern below: + +```cpp +template auto getQuadElement(Omega_h::Mesh &mesh) { + static_assert(ShapeOrder == 1); // extend as higher orders are added + if constexpr (ShapeOrder == 1) { + struct result { + MeshField::LinearQuadShape shp; + LinearQuadToVertexField map; + }; + return result{MeshField::LinearQuadShape(), + LinearQuadToVertexField(mesh)}; + } +} +``` + +Callers obtain the shape and mapping via structured bindings: + +```cpp +const auto [shp, map] = MeshField::Omegah::getQuadElement<1>(mesh); +MeshField::FieldElement fes(mesh.nelems(), field, shp, map); +``` + +--- + +## Step 4 – Add Integration Support + +`src/MeshField_Integrate.hpp` provides quadrature rules through the +`getIntegration()` template function. To integrate over +your new topology you must: + +### 4a – Define an `EntityIntegration` Class + +Derive from `EntityIntegration` where `D = meshEntDim` of your shape. +Implement at least one `Integration` inner class that returns a vector +of `IntegrationPoint` objects. Each point carries: + +| Field | Meaning | +|-------|---------| +| `param` | Parametric coordinates (reduced barycentric, length `D`) | +| `weight` | Quadrature weight | +| `dim` | Topological dimension of the entity the point is classified on | +| `idx` | Local entity index within the element | + +For the triangle, quadrature weights include the reference area factor +(\f$1/2\f$); for the tetrahedron they include \f$1/6\f$. Match the +convention used by `TriangleIntegration` and `TetrahedronIntegration` +already in the file. + +```cpp +class QuadIntegration : public EntityIntegration<2> { +public: + class N1 : public Integration<2> { + public: + int countPoints() const override { return 1; } + std::vector> getPoints() const override { + // 1-point Gauss rule for the reference quad [-1,1]^2, weight=4*(1/4)=1 + return {IntegrationPoint(Vector2{0.0, 0.0}, 1.0, 2, 0)}; + } + int getAccuracy() const override { return 1; } + }; + int countIntegrations() const override { return 1; } + Integration<2> const *getIntegration(int i) const override { + static N1 i1; + static Integration<2> *integrations[1] = {&i1}; + return integrations[i]; + } +}; +``` + +### 4b – Extend `getIntegration()` + +Add a branch to the existing `getIntegration` function: + +```cpp +template auto const getIntegration() { + if constexpr (topo == Triangle) { + return std::make_shared(); + } else if constexpr (topo == Tetrahedron) { + return std::make_shared(); + } else if constexpr (topo == Quad) { // <-- add this + return std::make_shared(); + } + fail("getIntegration does not support given topology\n"); +} +``` + +--- + +## Step 5 – Implement an Integrator (Optional) + +Derive from `MeshField::Integrator` and override `atPoints` to perform +the actual integration. `Integrator::process(FieldElement)` handles +quadrature point setup and calls `atPoints` with: + +| Argument | Type | Description | +|----------|------|-------------| +| `p` | `Kokkos::View` | Local parametric coordinates, shape `(numElems * numPts, dim)` | +| `w` | `Kokkos::View` | Quadrature weights, length `numElems * numPts` | +| `dV` | `Kokkos::View` | Jacobian determinants, length `numElems * numPts` | + +```cpp +template +class MyIntegrator : public MeshField::Integrator { + FieldElement &fes; + MeshField::Real result = 0.0; +public: + MyIntegrator(FieldElement &fe) : Integrator(/*order=*/1), fes(fe) {} + + void atPoints(Kokkos::View p, + Kokkos::View w, + Kokkos::View dV) override { + MeshField::Real local = 0.0; + Kokkos::parallel_reduce( + p.extent(0), + KOKKOS_LAMBDA(int i, MeshField::Real &sum) { sum += w(i) * dV(i); }, + local); + result += local; + } + + MeshField::Real getResult() const { return result; } +}; +``` + +Invoke it: + +```cpp +const auto [shp, map] = MeshField::Omegah::getQuadElement<1>(mesh); +MeshField::FieldElement fes(mesh.nelems(), coordField, shp, map); +MyIntegrator integrator(fes); +integrator.process(fes); +``` + +See `test/testCountIntegrator.cpp` for a working end-to-end example using +the existing triangle and tetrahedron shapes. + +--- + +## Checklist + +- [ ] Shape struct added to `src/MeshField_Shape.hpp` with all required members +- [ ] `getValues` and `getLocalGradients` validated at each node's parametric coords (partition-of-unity and Kronecker-delta checks) +- [ ] Omega\_h mapping struct added to `src/MeshField.hpp`, `namespace MeshField::Omegah` +- [ ] Vertex/edge ordering offset determined and verified against Omega\_h simplex templates +- [ ] Factory function `getElement` added or extended in `src/MeshField.hpp` +- [ ] `EntityIntegration` class added to `src/MeshField_Integrate.hpp` +- [ ] `getIntegration()` extended for the new topology +- [ ] Test added that constructs a mesh, builds a `FieldElement`, and calls `Integrator::process` From 133bcb6a06ecb40800c8d9cd1f06750458758216 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Thu, 6 Aug 2026 15:02:48 -0400 Subject: [PATCH 02/22] add to internal --- Doxyfile_internal.in | 1 + 1 file changed, 1 insertion(+) diff --git a/Doxyfile_internal.in b/Doxyfile_internal.in index 5fc54bbf..176665ed 100644 --- a/Doxyfile_internal.in +++ b/Doxyfile_internal.in @@ -919,6 +919,7 @@ WARN_LOGFILE = INPUT = \ @CMAKE_CURRENT_SOURCE_DIR@/docs/MeshFields.dox \ +@CMAKE_CURRENT_SOURCE_DIR@/docs/addingShapeFunctionsForOmegah.md \ @CMAKE_CURRENT_SOURCE_DIR@/src/CabanaController.hpp \ @CMAKE_CURRENT_SOURCE_DIR@/src/KokkosController.hpp \ @CMAKE_CURRENT_SOURCE_DIR@/src/MeshField_Defines.hpp \ From 492ae3b601944bc7224f3b227d5737f48f485b7b Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 7 Aug 2026 10:11:51 -0400 Subject: [PATCH 03/22] use snippet to reference code, don't hardcode examples --- Doxyfile.in | 2 +- Doxyfile_internal.in | 2 +- docs/addingShapeFunctionsForOmegah.md | 46 ++------------------------- src/MeshField_Shape.hpp | 2 ++ 4 files changed, 6 insertions(+), 46 deletions(-) diff --git a/Doxyfile.in b/Doxyfile.in index d6f1bce2..5dfa3e42 100644 --- a/Doxyfile.in +++ b/Doxyfile.in @@ -1057,7 +1057,7 @@ EXCLUDE_SYMBOLS = # that contain example code fragments that are included (see the \include # command). -EXAMPLE_PATH = +EXAMPLE_PATH = @CMAKE_CURRENT_SOURCE_DIR@ # If the value of the EXAMPLE_PATH tag contains directories, you can use the # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp and diff --git a/Doxyfile_internal.in b/Doxyfile_internal.in index 176665ed..2efdf40b 100644 --- a/Doxyfile_internal.in +++ b/Doxyfile_internal.in @@ -1065,7 +1065,7 @@ EXCLUDE_SYMBOLS = # that contain example code fragments that are included (see the \include # command). -EXAMPLE_PATH = +EXAMPLE_PATH = @CMAKE_CURRENT_SOURCE_DIR@ # If the value of the EXAMPLE_PATH tag contains directories, you can use the # EXAMPLE_PATTERNS tag to specify one or more wildcard pattern (like *.cpp and diff --git a/docs/addingShapeFunctionsForOmegah.md b/docs/addingShapeFunctionsForOmegah.md index 10ddcd04..790379bd 100644 --- a/docs/addingShapeFunctionsForOmegah.md +++ b/docs/addingShapeFunctionsForOmegah.md @@ -83,51 +83,9 @@ const Real L0 = 1 - xi[0] - xi[1]; assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol)); ``` -### Example Skeleton (Linear Quadrilateral) +### Example: `QuadraticTriangleShape` (from `src/MeshField_Shape.hpp`) -```cpp -struct LinearQuadShape { - static const size_t numNodes = 4; - static const size_t meshEntDim = 2; - constexpr static Mesh_Topology DofHolders[1] = {Vertex}; - constexpr static size_t Order = 1; - - KOKKOS_INLINE_FUNCTION - Kokkos::Array getNodeParametricCoords() const { - // clang-format off - return {-1,-1, // node 0 - 1,-1, // node 1 - 1, 1, // node 2 - -1, 1}; // node 3 - // clang-format on - } - - KOKKOS_INLINE_FUNCTION - Kokkos::Array getValues(Vector2 const &xi) const { - // standard bilinear shape functions - const Real xim = 1.0 - xi[0], xip = 1.0 + xi[0]; - const Real etam = 1.0 - xi[1], etap = 1.0 + xi[1]; - // clang-format off - return {0.25*xim*etam, - 0.25*xip*etam, - 0.25*xip*etap, - 0.25*xim*etap}; - // clang-format on - } - - KOKKOS_INLINE_FUNCTION - Kokkos::Array getLocalGradients(Vector2 const &xi) const { - const Real xim = 1.0 - xi[0], xip = 1.0 + xi[0]; - const Real etam = 1.0 - xi[1], etap = 1.0 + xi[1]; - // clang-format off - return {-0.25*etam, -0.25*xim, // dN0/dxi0, dN0/dxi1 - 0.25*etam, -0.25*xip, // dN1/dxi0, dN1/dxi1 - 0.25*etap, 0.25*xip, // dN2/dxi0, dN2/dxi1 - -0.25*etap, 0.25*xim}; // dN3/dxi0, dN3/dxi1 - // clang-format on - } -}; -``` +\snippet src/MeshField_Shape.hpp QuadraticTriangleShape --- diff --git a/src/MeshField_Shape.hpp b/src/MeshField_Shape.hpp index 8d5e9de4..e6a96432 100644 --- a/src/MeshField_Shape.hpp +++ b/src/MeshField_Shape.hpp @@ -269,6 +269,7 @@ struct LinearTriangleCoordinateShape { * Zienkiewicz, Taylor, and Zhu * 'The Finite Element Method: Its Basis and Fundamentals', 2013 */ +//! [QuadraticTriangleShape] struct QuadraticTriangleShape { static const size_t numNodes = 6; ///< Number of nodes (6 for quadratic triangle) static const size_t meshEntDim = 2; ///< Mesh entity dimension (2D) @@ -351,6 +352,7 @@ struct QuadraticTriangleShape { // clang-format on } }; +//! [QuadraticTriangleShape] /** * @brief Linear (P1) shape functions for 3D tetrahedral elements From f17256ee263dfac7b080c954cccae66b14f78c09 Mon Sep 17 00:00:00 2001 From: Cameron Smith Date: Fri, 7 Aug 2026 10:22:42 -0400 Subject: [PATCH 04/22] replace more hardcoded code --- docs/addingShapeFunctionsForOmegah.md | 130 +++----------------------- src/MeshField.hpp | 4 + src/MeshField_Integrate.hpp | 4 + test/testCountIntegrator.cpp | 4 + 4 files changed, 26 insertions(+), 116 deletions(-) diff --git a/docs/addingShapeFunctionsForOmegah.md b/docs/addingShapeFunctionsForOmegah.md index 790379bd..54a73813 100644 --- a/docs/addingShapeFunctionsForOmegah.md +++ b/docs/addingShapeFunctionsForOmegah.md @@ -140,44 +140,9 @@ struct. Validate with a unit test that evaluates the shape functions at the parametric coordinates of each node and checks that the value for that node is 1.0 and all others are 0.0. -### Example Skeleton (Linear Quad, vertices only) +### Example: `QuadraticTriangleToField` (from `src/MeshField.hpp`) -```cpp -struct LinearQuadToVertexField { - Omega_h::LOs quadVerts; - - LinearQuadToVertexField(Omega_h::Mesh &mesh) - : quadVerts(mesh.ask_elem_verts()) { - if (mesh.dim() != 2 || mesh.family() != OMEGA_H_HYPERCUBE) - MeshField::fail("Mesh must be 2D hypercube (quads)\n"); - } - - static constexpr KOKKOS_FUNCTION Kokkos::Array - getTopology() { - return {MeshField::Quad}; - } - - KOKKOS_FUNCTION MeshField::ElementToDofHolderMap - operator()(MeshField::LO nodeIdx, MeshField::LO compIdx, - MeshField::LO elem, MeshField::Mesh_Topology topo) const { - assert(topo == MeshField::Quad); - const auto quadDim = 2, vtxDim = 0; - // Determine the rotation offset by comparing Omega_h vertex ordering - // with your shape function's getNodeParametricCoords() ordering. - const auto localVtxIdx = - (Omega_h::simplex_down_template(quadDim, vtxDim, nodeIdx, -1) + OFFSET) % 4; - const auto stride = Omega_h::simplex_degree(quadDim, vtxDim); - const MeshField::LO vtx = quadVerts[elem * stride + localVtxIdx]; - return {0, compIdx, vtx, MeshField::Vertex}; - } -}; -``` - -Replace `OFFSET` with the value (0–3) that maps Omega\_h's vertex order -to the meshFields canonical order for your element. See the existing -`LinearTriangleToVertexField` (offset 2) and -`LinearTetrahedronToVertexField` (offset 3) in `src/MeshField.hpp` for -reference. +\snippet src/MeshField.hpp QuadraticTriangleToField --- @@ -185,26 +150,14 @@ reference. Add a new factory function (or extend an existing one) in `src/MeshField.hpp` inside `namespace MeshField::Omegah`. -The existing helpers follow the pattern below: +The existing `getTriangleElement` factory (from `src/MeshField.hpp`) shows the full pattern: -```cpp -template auto getQuadElement(Omega_h::Mesh &mesh) { - static_assert(ShapeOrder == 1); // extend as higher orders are added - if constexpr (ShapeOrder == 1) { - struct result { - MeshField::LinearQuadShape shp; - LinearQuadToVertexField map; - }; - return result{MeshField::LinearQuadShape(), - LinearQuadToVertexField(mesh)}; - } -} -``` +\snippet src/MeshField.hpp getTriangleElement Callers obtain the shape and mapping via structured bindings: ```cpp -const auto [shp, map] = MeshField::Omegah::getQuadElement<1>(mesh); +const auto [shp, map] = MeshField::Omegah::getTriangleElement<2>(mesh); MeshField::FieldElement fes(mesh.nelems(), field, shp, map); ``` @@ -234,43 +187,15 @@ For the triangle, quadrature weights include the reference area factor convention used by `TriangleIntegration` and `TetrahedronIntegration` already in the file. -```cpp -class QuadIntegration : public EntityIntegration<2> { -public: - class N1 : public Integration<2> { - public: - int countPoints() const override { return 1; } - std::vector> getPoints() const override { - // 1-point Gauss rule for the reference quad [-1,1]^2, weight=4*(1/4)=1 - return {IntegrationPoint(Vector2{0.0, 0.0}, 1.0, 2, 0)}; - } - int getAccuracy() const override { return 1; } - }; - int countIntegrations() const override { return 1; } - Integration<2> const *getIntegration(int i) const override { - static N1 i1; - static Integration<2> *integrations[1] = {&i1}; - return integrations[i]; - } -}; -``` +The existing `TriangleIntegration` (from `src/MeshField_Integrate.hpp`) shows the full pattern including multiple quadrature orders: + +\snippet src/MeshField_Integrate.hpp TriangleIntegration ### 4b – Extend `getIntegration()` -Add a branch to the existing `getIntegration` function: +Add a branch to the existing `getIntegration` function (current state in `src/MeshField_Integrate.hpp`): -```cpp -template auto const getIntegration() { - if constexpr (topo == Triangle) { - return std::make_shared(); - } else if constexpr (topo == Tetrahedron) { - return std::make_shared(); - } else if constexpr (topo == Quad) { // <-- add this - return std::make_shared(); - } - fail("getIntegration does not support given topology\n"); -} -``` +\snippet src/MeshField_Integrate.hpp getIntegration --- @@ -286,40 +211,13 @@ quadrature point setup and calls `atPoints` with: | `w` | `Kokkos::View` | Quadrature weights, length `numElems * numPts` | | `dV` | `Kokkos::View` | Jacobian determinants, length `numElems * numPts` | -```cpp -template -class MyIntegrator : public MeshField::Integrator { - FieldElement &fes; - MeshField::Real result = 0.0; -public: - MyIntegrator(FieldElement &fe) : Integrator(/*order=*/1), fes(fe) {} - - void atPoints(Kokkos::View p, - Kokkos::View w, - Kokkos::View dV) override { - MeshField::Real local = 0.0; - Kokkos::parallel_reduce( - p.extent(0), - KOKKOS_LAMBDA(int i, MeshField::Real &sum) { sum += w(i) * dV(i); }, - local); - result += local; - } - - MeshField::Real getResult() const { return result; } -}; -``` +`test/testCountIntegrator.cpp` provides a working example. The integrator class: -Invoke it: +\snippet test/testCountIntegrator.cpp CountIntegrator -```cpp -const auto [shp, map] = MeshField::Omegah::getQuadElement<1>(mesh); -MeshField::FieldElement fes(mesh.nelems(), coordField, shp, map); -MyIntegrator integrator(fes); -integrator.process(fes); -``` +And the function that wires together the element factory, `FieldElement`, integrator, and `process` call: -See `test/testCountIntegrator.cpp` for a working end-to-end example using -the existing triangle and tetrahedron shapes. +\snippet test/testCountIntegrator.cpp doRun --- diff --git a/src/MeshField.hpp b/src/MeshField.hpp index 7b24133d..dd541b86 100644 --- a/src/MeshField.hpp +++ b/src/MeshField.hpp @@ -124,6 +124,7 @@ struct LinearTetrahedronToVertexField { return {0, tetCompIdx, vtx, MeshField::Vertex}; // node, comp, ent, topo } }; +//! [QuadraticTriangleToField] struct QuadraticTriangleToField { Omega_h::LOs triVerts; Omega_h::LOs triEdges; @@ -185,6 +186,7 @@ struct QuadraticTriangleToField { return {0, triCompIdx, osh_ent, dofHolderTopo}; } }; +//! [QuadraticTriangleToField] struct QuadraticTetrahedronToField { Omega_h::LOs tetVerts; @@ -241,6 +243,7 @@ struct QuadraticTetrahedronToField { } }; +//! [getTriangleElement] template auto getTriangleElement(Omega_h::Mesh &mesh) { static_assert(ShapeOrder == 1 || ShapeOrder == 2); if constexpr (ShapeOrder == 1) { @@ -259,6 +262,7 @@ template auto getTriangleElement(Omega_h::Mesh &mesh) { QuadraticTriangleToField(mesh)}; } } +//! [getTriangleElement] template auto getTetrahedronElement(Omega_h::Mesh &mesh) { static_assert(ShapeOrder == 1 || ShapeOrder == 2); if constexpr (ShapeOrder == 1) { diff --git a/src/MeshField_Integrate.hpp b/src/MeshField_Integrate.hpp index 59ee529b..95ab549e 100644 --- a/src/MeshField_Integrate.hpp +++ b/src/MeshField_Integrate.hpp @@ -50,6 +50,7 @@ template class EntityIntegration { // points are given in the reduced parametric coordinates (xi0=L1, xi1=L2, // ...); the redundant barycentric coordinate L0 = 1-sum(xi) is omitted +//! [TriangleIntegration] class TriangleIntegration : public EntityIntegration<2> { public: class N1 : public Integration<2> { @@ -81,6 +82,7 @@ class TriangleIntegration : public EntityIntegration<2> { return integrations[i]; } }; +//! [TriangleIntegration] class TetrahedronIntegration : public EntityIntegration<3> { public: @@ -120,6 +122,7 @@ class TetrahedronIntegration : public EntityIntegration<3> { return integrations[i]; } }; +//! [getIntegration] template auto const getIntegration() { if constexpr (topo == Triangle) { return std::make_shared(); @@ -128,6 +131,7 @@ template auto const getIntegration() { } fail("getIntegration does not support given topology\n"); } +//! [getIntegration] template auto getIntegrationPoints(int order) { auto ip = getIntegration()->getAccurate(order)->getPoints(); return ip; diff --git a/test/testCountIntegrator.cpp b/test/testCountIntegrator.cpp index 763f61d2..781c179d 100644 --- a/test/testCountIntegrator.cpp +++ b/test/testCountIntegrator.cpp @@ -37,6 +37,7 @@ void setVertices(Omega_h::Mesh &mesh, AnalyticFunction func, ShapeField field) { setFieldAtVertices, "setFieldAtVertices"); } +//! [CountIntegrator] template class CountIntegrator : public MeshField::Integrator { private: @@ -56,7 +57,9 @@ class CountIntegrator : public MeshField::Integrator { count = fes.numMeshEnts; } }; +//! [CountIntegrator] +//! [doRun] template