Adding Shape Functions for Omega_h Element Topologies

Table of Contents

• Background
  • Terminology
  • Parametric Coordinates and Node Ordering
• Step 1 - Define the Shape Function Struct
  • Required Members
  • Checking Parametric Coordinates
  • Example: <tt>QuadraticTriangleShape</tt> (from <tt>src/MeshField_Shape.hpp</tt>)
• Step 2 - Define the Omega_h Node Mapping Struct
  • Vertex-Ordering Correction
  • Interface Requirements
  • Omega_h Connectivity APIs
  • Example: <tt>QuadraticTriangleToField</tt> (from <tt>src/MeshField.hpp</tt>)
• Step 3 - Register in the Element Factory
• Step 4 - Add an Accessor
• Step 5 - Extend <tt>CreateLagrangeField</tt>

Shape functions for an element topology are defined within a shape function struct (e.g., MeshField::LinearTriangleShape). Mapping between the ordering used for those shape functions and those used in the mesh database (currently just Omega_h) is through a mapping struct (see Step 2 - Define the Omega_h Node Mapping Struct below). DOF values are stored in MeshField::Field which needs to be associated with a shape function and mapping via a ShapeField to provide the needed per entity/dof access operators.

Note, as of MeshFields version 1.1.0 only interpolating Lagrange linear and quadratic triangles and tetrahedra are supported. As such, there may be deficiencies in the interface to support other types of shape functions. Please use https://github.com/SCOREC/meshFields/issues to report any problems or ask questions.

This guide walks through every step required to add a new shape function for an Omega_h element topology (edge, triangle, quadrilateral, tetrahedron, etc.). The steps are:

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 an Accessor (if needed) in src/MeshField_ShapeField.hpp
5. Extend **CreateLagrangeField** in src/MeshField_ShapeField.hpp

---

Background

Terminology

Refer to the Nomenclature section on the main page.

Parametric Coordinates and Node Ordering

MeshFields defines shape functions in a canonical parametric coordinate system that is independent of Omega_h. The node mapping struct (Step 2) bridges the two conventions.

The MeshFields canonical parametric coordinates and node ordering for linear and quadratic triangles and tetrahedrons follows, "The Finite Element Method: Its Basis and Fundamentals", 2013, Zienkiewicz, Taylor, and Zhu.

MeshFields uses parametric coordinates to specify a location within an element of dimension . The redundant coordinate is omitted:

Topology	Parametric coords	Range
Edge (1D)
Triangle (2D)		,
Tet (3D)		,

---

Step 1 - Define the Shape Function Struct

Back to top

Add a new struct to src/MeshField_Shape.hpp inside namespace MeshField.

Required Members

Member	Description
numNodes	Total nodes per element
meshEntDim	Parametric space dimension
Order	Polynomial order
DofHolders	Entity types that hold DOFs
getNodeParametricCoords()	Returns a flat array of node coordinates (length numNodes * meshEntDim)
getValues(xi)	Returns an array of shape function values at xi (length numNodes)
getLocalGradients(xi)	Returns a flat array of gradients at xi (length meshEntDim * numNodes, row-major: )

Optional (needed by quadratic and higher-order shapes):

Member	Description
NumDofHolders	Count of entities of each type in DofHolders
DofsPerHolder	DOFs per entity for each type in DofHolders

Checking Parametric Coordinates

Use the helper functions declared at the top of MeshField_Shape.hpp to assert that incoming parametric coordinates are in range. This is only a check for outliers as MeshField::ParametricCoordTol is a loose tolerance.

assert(eachGreaterThanOrEqual(xi, 0.0, ParametricCoordTol));
assert(eachLessThanOrEqual(xi, 1.0, ParametricCoordTol));
const Real L0 = 1 - xi[0] - xi[1];
assert(greaterThanOrEqual(L0, 0.0, ParametricCoordTol));

Example: <tt>QuadraticTriangleShape</tt> (from <tt>src/MeshField_Shape.hpp</tt>)

struct QuadraticTriangleShape {
  static const size_t numNodes = 6;                              
  static const size_t meshEntDim = 2;                            
  constexpr static Mesh_Topology DofHolders[2] = {Vertex, Edge}; 
  constexpr static size_t NumDofHolders[2] = {3, 3};             
  constexpr static size_t DofsPerHolder[2] = {1, 1};             
  constexpr static size_t Order = 2;                             
 
  KOKKOS_INLINE_FUNCTION
  Kokkos::Array<Real, numNodes * meshEntDim> getNodeParametricCoords() const {
    // clang-format off
    return {
      //nodes at vertices
      0   , 0   , //node 0
      1   , 0   , //node 1
      0   , 1   , //...
      //nodes at middle of edges
      0.5 , 0   ,
      0.5 , 0.5 ,
      0   , 0.5   //node 5
    };
    // clang-format on
  }
 
  KOKKOS_INLINE_FUNCTION
  Kokkos::Array<Real, numNodes> getValues(Vector2 const &xi) const {
    assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol));
    assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol));
    const Real L0 = 1 - xi[0] - xi[1];
    assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol));
    assert(lessThanOrEqual(L0,1.0,ParametricCoordTol));
    const Real L1 = xi[0];
    const Real L2 = xi[1];
    // clang-format off
    return {L0 * (2 * L0 - 1),
            L1 * (2 * L1 - 1),
            L2 * (2 * L2 - 1),
            4 * L1 * L0,
            4 * L1 * L2,
            4 * L2 * L0};
    // clang-format on
  }
 
  KOKKOS_INLINE_FUNCTION
  Kokkos::Array<Real, meshEntDim * numNodes>
  getLocalGradients(Vector2 const &xi) const {
    assert(eachLessThanOrEqual(xi,1.0,ParametricCoordTol));
    assert(eachGreaterThanOrEqual(xi,0.0,ParametricCoordTol));
    const Real L0 = 1 - xi[0] - xi[1];
    assert(greaterThanOrEqual(L0,0.0,ParametricCoordTol));
    assert(lessThanOrEqual(L0,1.0,ParametricCoordTol));
    const Real L1 = xi[0];
    const Real L2 = xi[1];
    // clang-format off
    return {
      -4*L0+1   , -4*L0+1   ,
      4*L1-1    , 0         ,
      0         , 4*L2-1    ,
      4*(L0-L1) , -4*L1     ,
      4*L2      , 4*L1      ,
      -4*L2     , 4*(L0-L2)
    };
    // clang-format on
  }
};

---

Step 2 - Define the Omega_h Node Mapping Struct

Back to top

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 on-process mesh entity indices (as numbered by Omega_h) that hold the corresponding DOFs.

Vertex-Ordering Correction

As previously stated, Omega_h numbers vertices and edges within an element differently from the MeshFields canonical ordering used by the shape functions. The Omega_h canonical orderings for simplices are depicted below:

Omega_h triangle ordering.

Omega_h tetrahedron ordering.

Meshes composed of hypercubes (quadrilaterals and hexahedrons), pyramids, and prisms are supported by Omega_h, but they cannot be adapted. Version 1.1.0 of MeshFields only supports Omega_h simplices.

The existing linear-triangle and linear-tetrahedron mappings correct for the difference in Omega_h vs MeshFields ordering with a cyclic rotation.

// 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 and DOF holders by comparing the Omega_h canonical vertex/edge ordering (see figures above based on Omega_h_simplex.hpp and the simplex_down_template function) with the node ordering established by getNodeParametricCoords() in your shape struct. It is a good idea to add 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.

Interface Requirements

Member	Description
Constructor (Omega_h::Mesh &)	Cache connectivity arrays from Omega_h; validate mesh family/dimension
getTopology()	Return an array of mesh topologies (vertex, edge, triangle, etc...) this mapping applies to
operator()(LO nodeIndex_in, LO componentIndex_in, LO elementIndex_in, Mesh_Topology elementTopo_in)	Return mapping tuple - see MeshField::ElementToDofHolderMap

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

Example: <tt>QuadraticTriangleToField</tt> (from <tt>src/MeshField.hpp</tt>)

struct QuadraticTriangleToField {
  Omega_h::LOs triVerts;
  Omega_h::LOs triEdges;
  QuadraticTriangleToField(Omega_h::Mesh &mesh)
      : triVerts(mesh.ask_elem_verts()),
        triEdges(mesh.ask_down(mesh.dim(), 1).ab2b) {
    if (mesh.dim() != 2 && mesh.family() != OMEGA_H_SIMPLEX) {
      MeshField::fail(
          "The mesh passed to %s must be 2D and simplex (triangles)\n",
          __func__);
    }
  }
 
  static constexpr KOKKOS_FUNCTION Kokkos::Array<MeshField::Mesh_Topology, 1>
  getTopology() {
    return {MeshField::Triangle};
  }
 
  KOKKOS_FUNCTION MeshField::ElementToDofHolderMap
  operator()(MeshField::LO triNodeIdx, MeshField::LO triCompIdx,
             MeshField::LO tri, MeshField::Mesh_Topology topo) const {
    assert(topo == MeshField::Triangle);
    // Omega_h has no concept of nodes so we can define the map from
    // triNodeIdx to the dof holder index
    const MeshField::LO triNode2DofHolder[6] = {
        /*vertices*/ 0, 1, 2,
        /*edges*/ 0,    1, 2};
    const MeshField::Mesh_Topology triNode2DofHolderTopo[6] = {
        /*vertices*/
        MeshField::Vertex, MeshField::Vertex, MeshField::Vertex,
        /*edges*/
        MeshField::Edge, MeshField::Edge, MeshField::Edge};
    const auto dofHolderIdx = triNode2DofHolder[triNodeIdx];
    const auto dofHolderTopo = triNode2DofHolderTopo[triNodeIdx];
    // Given the topo index and type find the Omega_h vertex or edge index that
    // bounds the triangle
    Omega_h::LO osh_ent;
    if (dofHolderTopo == MeshField::Vertex) {
      const auto triDim = 2;
      const auto vtxDim = 0;
      const auto ignored = -1;
      const auto localVtxIdx = (Omega_h::simplex_down_template(
                                    triDim, vtxDim, dofHolderIdx, ignored) +
                                2) %
                               3;
      const auto triToVtxDegree = Omega_h::simplex_degree(triDim, vtxDim);
      osh_ent = triVerts[(tri * triToVtxDegree) + localVtxIdx];
    } else if (dofHolderTopo == MeshField::Edge) {
      const auto triDim = 2;
      const auto edgeDim = 1;
      const auto triToEdgeDegree = Omega_h::simplex_degree(triDim, edgeDim);
      // passing dofHolderIdx as Omega_h_simplex.hpp does not provide
      // a function that maps a triangle and edge index to a 'canonical' edge
      // index. This may need to be revisited...
      osh_ent = triEdges[(tri * triToEdgeDegree) + (dofHolderIdx + 2) % 3];
    } else {
      assert(false);
    }
    return {0, triCompIdx, osh_ent, dofHolderTopo};
  }
};

---

Step 3 - Register in the Element Factory

Back to top

Add a new factory function (or extend an existing one) in src/MeshField.hpp inside namespace MeshField::Omegah. The existing getTriangleElement factory (from src/MeshField.hpp) shows the full pattern:

template <int ShapeOrder> auto getTriangleElement(Omega_h::Mesh &mesh) {
  static_assert(ShapeOrder == 1 || ShapeOrder == 2);
  if constexpr (ShapeOrder == 1) {
    struct result {
      MeshField::LinearTriangleShape shp;
      LinearTriangleToVertexField map;
    };
    return result{MeshField::LinearTriangleShape(),
                  LinearTriangleToVertexField(mesh)};
  } else if constexpr (ShapeOrder == 2) {
    struct result {
      MeshField::QuadraticTriangleShape shp;
      QuadraticTriangleToField map;
    };
    return result{MeshField::QuadraticTriangleShape(),
                  QuadraticTriangleToField(mesh)};
  }
}

Callers obtain the shape and mapping via structured bindings:

const auto [shp, map] = MeshField::Omegah::getTriangleElement<2>(mesh);
MeshField::FieldElement fes(mesh.nelems(), field, shp, map);

---

Step 4 - Add an Accessor

Back to top

An Accessor is a templated struct that provides operator()(entity, node, component, topology) for reading and writing DOF values from the underlying field storage. It is passed as a variadic template argument (Mixins...) to ShapeField.

Two built-in Accessors cover the common cases:

LinearAccessor – for shapes whose DOFs live only at vertices. All calls to operator() are forwarded to a single vertex field regardless of topology:

template <typename VtxAccessor> struct LinearAccessor {
  constexpr static const Mesh_Topology topo[1] = {Vertex};
  VtxAccessor vtxField;
  using BaseType = typename VtxAccessor::BaseType;
 
  KOKKOS_FUNCTION
  auto &operator()(int entity, int node, int component, Mesh_Topology t) const {
    if (t != Vertex) {
      Kokkos::printf("%d is not a support topology\n", t);
      assert(false);
    }
    return vtxField(entity, node, component);
  }
};

QuadraticAccessor – for shapes whose DOFs live at both vertices and edges. The operator() dispatches to the vertex field or edge field based on the topology argument:

template <typename VtxAccessor, typename EdgeAccessor>
struct QuadraticAccessor {
  constexpr static const Mesh_Topology topo[2] = {Vertex, Edge};
  VtxAccessor vtxField;
  EdgeAccessor edgeField;
  using BaseType = typename VtxAccessor::BaseType;
 
  KOKKOS_FUNCTION
  auto &operator()(int entity, int node, int component, Mesh_Topology t) const {
    if (t != Vertex && t != Edge) {
      Kokkos::printf("%d is not a support topology\n", t);
      assert(false);
    }
    return (t == Vertex) ? vtxField(entity, node, component)
                         : edgeField(entity, node, component);
  }
};

If your shape introduces DOFs at a new entity type (e.g., triangle faces for a cubic shape), define a new Accessor following the same pattern:

1. Add a field member (as done in QuadraticAccessor for EdgeAccessor edgeField;,) for each new entity type.
2. Extend operator() with a branch for the added topology.
3. List all supported topologies in the topo array.

---

Step 5 - Extend <tt>CreateLagrangeField</tt>

Back to top

CreateLagrangeField<ExecutionSpace, Controller, DataType, order, dim, numComp>(meshInfo) (in src/MeshField_ShapeField.hpp) is the user-facing factory that allocates field storage and assembles a ShapeField. Add a new if constexpr branch for your shape's order and dim. The steps within that branch are:

1. Validate meshInfo counts for every entity type that holds DOFs (e.g., meshInfo.numVtx, meshInfo.numEdge). Call fail() if a required count is zero.
2. Allocate storage – construct a Controller sized for each DOF-holding entity type. Each field in the controller needs entries for numEntities * dofsPerHolder * numComp.
3. Create fields – call MeshField::makeField<Ctrlr, N>(ctrl) for each field index N to obtain typed slice objects.
4. Construct the Accessor – pass the slice objects to your Accessor's constructor (e.g., QuadraticAccessor{vtxField, edgeField}).
5. Return FieldWithController<Ctrlr, ShapeField<numComp, YourShape, YourAccessor>>{ctrl, ShapeField(meshInfo, accessor)}.

See the existing order == 2 && dim == 2 branch (quadratic triangle) as a concrete reference:

    if (meshInfo.numVtx <= 0) {
      fail("mesh has no vertices\n");
    }
    if (meshInfo.numEdge <= 0) {
      fail("mesh has no edges\n");
    }
#ifdef MESHFIELDS_ENABLE_CABANA
    using Ctrlr = std::conditional_t<
        std::is_same_v<
            Controller<ExecutionSpace, MemorySpace, DataType>,
            MeshField::CabanaController<ExecutionSpace, MemorySpace, DataType>>,
        Controller<ExecutionSpace, MemorySpace, DataType[1][numComp],
                   DataType[1][numComp]>,
        Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>>;
    // 1 dof with 1 comp per vtx/edge
    auto createController = [](auto numVtx, auto numEdge) {
      if constexpr (std::is_same_v<
                        Controller<ExecutionSpace, MemorySpace, DataType>,
                        MeshField::CabanaController<ExecutionSpace, MemorySpace,
                                                    DataType>>) {
        return Ctrlr({numVtx, numEdge});
      } else {
        return Ctrlr({/*field 0*/ numVtx, 1, numComp,
                      /*field 1*/ numEdge, 1, numComp});
      }
    };
    Ctrlr kk_ctrl = createController(meshInfo.numVtx, meshInfo.numEdge);
#else
    using Ctrlr =
        Controller<MemorySpace, ExecutionSpace, DataType ***, DataType ***>;
    Ctrlr kk_ctrl({/*field 0*/ meshInfo.numVtx, 1, numComp,
                   /*field 1*/ meshInfo.numEdge, 1, numComp});
#endif
    auto vtxField = MeshField::makeField<Ctrlr, 0>(kk_ctrl);
    auto edgeField = MeshField::makeField<Ctrlr, 1>(kk_ctrl);
    using QA = QuadraticAccessor<decltype(vtxField), decltype(edgeField)>;
    // clang-format off
    using QuadraticLagrangeShapeField = std::conditional_t<
        dim == 3,
        ShapeField<numComp, QuadraticTetrahedronShape, QA>,
        ShapeField<numComp, QuadraticTriangleShape, QA>>;
    // clang-format on
    QuadraticLagrangeShapeField qlsf(meshInfo, {vtxField, edgeField});
    return FieldWithController<Ctrlr, QuadraticLagrangeShapeField>{kk_ctrl,
                                                                    qlsf};