MeshFields
GPU accelerated mesh-based fields
Adding Shape Functions for Omega_h Element Topologies

Table of Contents

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 $d-1$ parametric coordinates to specify a location within an element of dimension $d$. The redundant coordinate $L_0 = 1 - \sum \xi_i$ is omitted:

Topology Parametric coords Range
Edge (1D) $\xi$ $[-1, 1]$
Triangle (2D) $(\xi_0,\xi_1)$ $[0,1]^2$, $\xi_0+\xi_1 \le 1$
Tet (3D) $(\xi_0,\xi_1,\xi_2)$ $[0,1]^3$, $\sum \xi_i \le 1$

Diagrams of the vertex and edge node ordering, with parametric coordinates, for the linear and quadratic triangle and tetrahedron are in docs/Ordering.tex.

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: $[\partial N_0/\partial\xi_0, \partial N_0/\partial\xi_1, \ldots, \partial N_d/\partial\xi_0, \partial N_d/\partial\xi_1]$)

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.

Mapping Omega_h Vertices and Edges to Nodes

The Omega_h canonical orderings for simplices are depicted below:

Omega_h triangle ordering

Omega_h triangle ordering.

Omega_h tetrahedron 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.

Vertices. The MeshFields vertex node ordering for triangles and tetrahedra matches the Omega_h element-to-vertex ordering, so vertex node $i$ maps directly to local vertex $i$ with no rotation:

// For triangles (triDim=2) and tetrahedra (tetDim=3), vtxDim=0:
const auto localVtxIdx =
Omega_h::simplex_down_template(elemDim, vtxDim, nodeIdx, /*ignored=*/-1);

Edges. Edge nodes do not necessarily follow the Omega_h element-to-edge ordering, so the mapping struct uses a lookup table from node index to the Omega_h local edge index. For the quadratic triangle the MeshFields edges (E0: V0-V1, E1: V1-V2, E2: V0-V2) coincide with the Omega_h edges, so the table is the identity. For the quadratic tetrahedron the MeshFields edges (nodes 4-9: V0-V1, V0-V2, V0-V3, V1-V2, V2-V3, V1-V3) are reordered to the Omega_h edges (0: V0-V1, 1: V1-V2, 2: V2-V0, 3: V0-V3, 4: V1-V3, 5: V2-V3):

// from QuadraticTetrahedronToField
const MeshField::LO tetNode2DofHolder[10] = {0, 1, 2, 3, 0, 2, 3, 1, 5, 4};

You must determine the correct vertex and edge mapping 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.

Jacobian Orientation

FieldElement returns signed Jacobian determinants for 2D and 3D elements. The integrator (src/MeshField_Integrate.hpp) uses the determinants as-is, and debug builds assert that each one is non-negative. The vertex ordering produced by your mapping must therefore follow the right-hand rule so that the Jacobian determinant of a valid (un-inverted) Omega_h element is positive. An incorrect vertex permutation results in a negative determinant and a failed assertion in debug builds.

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);
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];
} else {
assert(false);
}
return {0, triCompIdx, osh_ent, dofHolderTopo};
}
};
MeshField::ElementToDofHolderMap
Supports mapping between mesh (i.e., Omega_h) ordering and MeshFields ordering.
Definition: MeshField_Element.hpp:101

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)};
}
}
MeshField::LinearTriangleShape
Linear (P1) shape functions for 2D triangular elements.
Definition: MeshField_Shape.hpp:163
MeshField::QuadraticTriangleShape
Quadratic (P2) shape functions for 2D triangular elements.
Definition: MeshField_Shape.hpp:285

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);
MeshField::FieldElement
Supports the evaluation of a field, and other per-element operations, given the definition of the ele...
Definition: MeshField_Element.hpp:149

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, Ctrlr, 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, Ctrlr, QA>,
ShapeField<numComp, QuadraticTriangleShape, Ctrlr, QA>>;
// clang-format on
QuadraticLagrangeShapeField qlsf(meshInfo, {vtxField, edgeField});
return FieldWithController<Ctrlr, QuadraticLagrangeShapeField>{kk_ctrl,
qlsf};
Generated by doxygen 1.9.1