tensormesh.element

The tensormesh.element module is organised in three tiers, each aimed at a different audience.

1. Transformation — the only thing most users call directly. If you are writing a new assembler or computing shape values / gradients / J×W at quadrature points, this is your interface.

2. Reference-shape catalog — the seven reference shapes and the string-↔-class registry. Read this when you need to know which type string corresponds to which class, what its facet shape is, or how its node ordering compares to Gmsh / VTK.

3. Extension internals — basis builders, quadrature factories, facet-normal helpers, and the Polynomial machinery. Read this only if you are adding a new element type or quadrature rule.

Transformation

The workhorse: given a mesh’s points + connectivity + element type, a Transformation caches the per-element Jacobian, shape-function values, shape-function gradients, quadrature points/weights, and their facet analogues. Every built-in assembler in tensormesh.assemble ultimately reads from a Transformation; custom assemblers should too.

class Transformation(points: Tensor, elements: Tensor, element_type: str, quadrature_order: int = 2)

Bases: Module

Reference-to-physical mapping for one element-type block of a mesh.

Holds the per-element Jacobian, shape-function values, shape-function gradients, quadrature points/weights, and their facet analogues — all the geometric ingredients an assembler needs to integrate a weak form. Subclass of torch.nn.Module, so every cached tensor is a buffer and follows the module’s .to()/.cuda()/.double().

Examples

import tensormesh as tm

mesh = tm.Mesh.gen_rectangle(chara_length=0.1)        # triangle mesh
transform = tm.Transformation(
    mesh.points,
    mesh.cells["triangle"],
    element_type="triangle",
    quadrature_order=2,
)

# Geometric quantities at quadrature points
J   = transform.jacobian      # [n_elem, n_q, dim, dim]
phi = transform.shape_val     # [n_q, n_basis]
gphi= transform.shape_grad    # [n_elem, n_q, n_basis, dim]
jxw = transform.JxW           # [n_elem, n_q]  (|det J| · w)

elements

Element connectivity tensor of shape \([N_e, N_b]\) where:

• \(N_e\) = number of elements

• \(N_b\) = number of basis functions per element

Type:

Tensor

points

Point coordinates tensor of shape \([N_p, D]\) where:

• \(N_p\) = number of points

• \(D\) = spatial dimension

Type:

Tensor

element

Element type class (e.g. Triangle, Tetrahedron)

Type:

Type[Element]

basis_order

Order of basis functions

Type:

int

quadrature_order

Order of quadrature rule

Type:

int

__init__(points: Tensor, elements: Tensor, element_type: str, quadrature_order: int = 2)

Initialize transformation.

Examples

# Single linear triangle
points   = torch.tensor([[0., 0.], [1., 0.], [0., 1.]])
elements = torch.tensor([[0, 1, 2]])
transform = Transformation(points, elements, "triangle", quadrature_order=2)

Parameters:

• points (Tensor) – Point coordinates tensor of shape \([N_p, D]\).

• elements (Tensor) – Element connectivity tensor of shape \([N_e, N_b]\).

• element_type (str) – Element type string (e.g. "triangle", "triangle6", "tetra", "hexahedron27"). See tensormesh.element_types for the full list.

• quadrature_order (int, optional) – Order of the quadrature rule, by default 2.

update_points(points: Tensor)

Replace the cached node coordinates and reset derived buffers.

Use this after deforming the mesh — e.g. in a moving-mesh or topology-optimization loop — when the connectivity does not change but the geometry does.

Examples

new_points = transform.points + displacement
transform.update_points(new_points)

Parameters:

points (Tensor) – Point coordinates tensor of shape \([N_p, D]\) where \(N_p\) is the number of points and \(D\) the spatial dimension.

property dim: int

Get the spatial dimension of the mesh.

Returns:

The spatial dimension \(D\) (1, 2, or 3)

Return type:

int

property device: device

Get the device of the mesh tensors.

Returns:

The device (CPU/GPU) where the mesh tensors are stored

Return type:

device

property dtype: dtype

Get the data type of the mesh tensors.

Returns:

The data type (e.g. torch.float32) of the mesh tensors

Return type:

dtype

property element_coords: Tensor

Get the physical coordinates of element nodes.

Examples

transform = Transformation(points, elements, "triangle")
coords = transform.element_coords
# coords[e]    -> [n_basis, dim] coordinates of all nodes in element e
# coords[e, b] -> [dim] coordinate of basis node b in element e

Returns:

Element node coordinates tensor of shape \([N_e, N_b, D]\) where:

• \(N_e\) = number of elements

• \(N_b\) = number of basis functions per element

• \(D\) = spatial dimension

Return type:

Tensor

property n_elements: int

Get the number of elements in the mesh.

Returns:

Number of elements \(N_e\)

Return type:

int

property basis: Tensor

Get the basis function nodes in reference coordinates.

Examples

transform = Transformation(points, elements, "triangle", basis_order=2)
basis = transform.basis  # Get quadratic basis nodes
# For triangles: 6 nodes (3 vertices + 3 edge midpoints)

Returns:

Basis node coordinates tensor of shape \([N_b, D]\) where:

• \(N_b\) = number of basis functions

• \(D\) = spatial dimension

Return type:

Tensor

property n_basis: int

Get the number of basis functions per element.

Returns:

Number of basis functions \(N_b\)

Return type:

int

property quadrature: Tuple[Tensor, Tensor]

Get quadrature weights and points.

Examples

transform = Transformation(points, elements, "triangle", quadrature_order=2)
weights, points = transform.quadrature
# weights: quadrature weights for numerical integration
# points: quadrature point coordinates in reference element

Returns:

• Quadrature weights tensor of shape \([N_q]\)

• Quadrature points tensor of shape \([N_q, D]\)

where:

• \(N_q\) = number of quadrature points

• \(D\) = spatial dimension

Return type:

Tuple[Tensor, Tensor]

property n_quadrature: int

Get the number of quadrature points per element.

Returns:

Number of quadrature points \(N_q\)

Return type:

int

property shape_val: Tensor

Get shape function values at quadrature points.

The shape functions \(N_i(\boldsymbol{\xi})\) are evaluated at the quadrature points in the reference element.

Examples

transform = Transformation(points, elements, "triangle")
N = transform.shape_val
# N[i,j] is value of jth shape function at ith quadrature point

Returns:

Shape function values tensor of shape \([N_q, N_b]\) where:

• \(N_q\) = number of quadrature points

• \(N_b\) = number of basis functions

Return type:

Tensor

property shape_grad: Tensor

Get shape function gradients at quadrature points.

The shape function gradients \(\nabla N_i(\mathbf{x})\) are evaluated at the quadrature points in physical coordinates.

Examples

transform = Transformation(points, elements, "triangle")
dN = transform.shape_grad
# dN[e,q,i,d] is derivative of ith shape function w.r.t.
# dth coordinate at qth quadrature point in eth element

Returns:

Shape function gradients tensor of shape \([N_e, N_q, N_b, D]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

• \(N_b\) = number of basis functions

• \(D\) = spatial dimension

Return type:

Tensor

property jacobian: Tensor

Get Jacobian matrices at basis points,

The Jacobian matrix \(\mathbf{J} = \frac{\partial\mathbf{x}}{\partial\boldsymbol{\xi}}\) maps derivatives from reference to physical coordinates.

The Jacobian matrix maps derivatives between reference and physical coordinates:

\[\frac{\partial}{\partial x_i} = \sum_{j=1}^D J^{-1}_{ij} \frac{\partial}{\partial \xi_j}\]

where \(J_{ij} = \frac{\partial x_i}{\partial \xi_j}\) and \(D\) is the spatial dimension.

Examples

transform = Transformation(points, elements, "triangle")
J = transform.jacobian
# J[e,q,:,:] is Jacobian matrix at qth quadrature point
# in eth element

Returns:

Jacobian matrices tensor of shape \([N_e, N_q, D, D]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

• \(D\) = spatial dimension

Return type:

Tensor

property facets: Tensor | Tuple[Tensor, Tensor]

Get element facet connectivity.

For elements with uniform facets (e.g. all triangles or all quads), returns a single tensor. For elements with mixed facets (e.g. prism or pyramid), returns separate tensors for triangular and quadrilateral facets.

Examples

# For a tetrahedral mesh (uniform triangular facets)
transform = Transformation(points, elements, "tetra")
facets = transform.facets  # Single tensor of triangular facets

# For a prism mesh (mixed tri/quad facets)
transform = Transformation(points, elements, "wedge")
tri_facets, quad_facets = transform.facets  # Separate tensors

Returns:

• If element has uniform facets –

  torch.Tensor

  Facet connectivity tensor of shape \([N_f, N_b]\) where:

  • \(N_f\) = number of facets

  • \(N_b\) = number of basis functions per facet

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Triangular facet tensor of shape \([N_{tf}, N_{tb}]\)

  • Quadrilateral facet tensor of shape \([N_{qf}, N_{qb}]\)

  where:

  • \(N_{tf}\) = number of triangular facets

  • \(N_{tb}\) = number of basis functions per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets

  • \(N_{qb}\) = number of basis functions per quadrilateral facet

property facet_quadrature: Tuple[Tensor, Tensor] | Tuple[Tensor, Tensor, Tensor, Tensor]

Get quadrature weights and points for element facets.

Provides quadrature rules for numerical integration over element facets. For mixed facet elements, returns separate quadrature rules for triangular and quadrilateral facets.

Examples

# For tetrahedral mesh (triangular facets)
transform = Transformation(points, elements, "tetra")
weights, points = transform.facet_quadrature

# For prism mesh (mixed facets)
transform = Transformation(points, elements, "wedge")
tri_w, tri_p, quad_w, quad_p = transform.facet_quadrature

Returns:

• If element has uniform facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Weights tensor of shape \([N_f, N_q]\)

  • Points tensor of shape \([N_f, N_q, D]\)

  where:

  • \(N_f\) = number of facets

  • \(N_q\) = number of quadrature points per facet

  • \(D\) = spatial dimension

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor, torch.Tensor, torch.Tensor]

  • Triangular facet weights of shape \([N_{tf}, N_{tq}]\)

  • Triangular facet points of shape \([N_{tf}, N_{tq}, D]\)

  • Quadrilateral facet weights of shape \([N_{qf}, N_{qq}]\)

  • Quadrilateral facet points of shape \([N_{qf}, N_{qq}, D]\)

  where:

  • \(N_{tf}\) = number of triangular facets

  • \(N_{tq}\) = quadrature points per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets

  • \(N_{qq}\) = quadrature points per quadrilateral facet

  • \(D\) = spatial dimension

property facet_shape_val: Tensor | Tuple[Tensor, Tensor]

Get shape function values at facet quadrature points.

Evaluates the element shape functions at the quadrature points on each facet. For mixed facet elements, returns separate evaluations for triangular and quadrilateral facets.

The shape functions \(\phi_i(\mathbf{x})\) are evaluated at quadrature points on each facet:

\[\phi_i(\mathbf{x}_q^f) = \phi_i(\mathbf{F}^f(\mathbf{r}_q))\]

where:

• \(\mathbf{x}_q^f\) = physical coordinates of quadrature point \(q\) on facet \(f\)

• \(\mathbf{F}^f\) = mapping from reference to physical facet

• \(\mathbf{r}_q\) = reference coordinates of quadrature point \(q\)

Examples

# For tetrahedral mesh (triangular facets)
transform = Transformation(points, elements, "tetra")
shape_vals = transform.facet_shape_val

# For prism mesh (mixed facets)
transform = Transformation(points, elements, "wedge")
tri_vals, quad_vals = transform.facet_shape_val

Returns:

• If element has uniform facets –

  torch.Tensor

  Shape function values tensor of shape \([N_f, N_q, N_b]\) where:

  • \(N_f\) = number of facets

  • \(N_q\) = number of quadrature points per facet

  • \(N_b\) = number of basis functions

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Triangular facet values of shape \([N_{tf}, N_{tq}, N_b]\)

  • Quadrilateral facet values of shape \([N_{qf}, N_{qq}, N_b]\)

  where:

  • \(N_{tf}\) = number of triangular facets

  • \(N_{tq}\) = quadrature points per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets

  • \(N_{qq}\) = quadrature points per quadrilateral facet

  • \(N_b\) = number of basis functions

property facet_shape_grad: Tensor | Tuple[Tensor, Tensor]

Get shape function gradients at facet quadrature points.

Evaluates the gradients of element shape functions at quadrature points on each facet. For mixed facet elements, returns separate gradient evaluations for triangular and quadrilateral facets.

The shape function gradients \(\nabla\phi_i\) are evaluated at quadrature points on each facet. For a facet \(f\), the gradients are computed as:

\[\nabla\phi_i(\mathbf{x}_q^f) = \mathbf{J}^{-T}(\mathbf{x}_q^f) \nabla\hat{\phi}_i(\boldsymbol{\xi}_q^f)\]

where:

• \(\mathbf{x}_q^f\) is the physical coordinate of quadrature point \(q\) on facet \(f\)

• \(\boldsymbol{\xi}_q^f\) is the reference coordinate of the quadrature point

• \(\mathbf{J}\) is the Jacobian matrix

• \(\hat{\phi}_i\) is the reference element shape function

Examples

# For tetrahedral mesh (triangular facets)
transform = Transformation(points, elements, "tetra")
grads = transform.facet_shape_grad

# For prism mesh (mixed facets)
transform = Transformation(points, elements, "wedge")
tri_grads, quad_grads = transform.facet_shape_grad

Returns:

• If element has uniform facets –

  torch.Tensor

  Shape function gradients tensor of shape \([N_f, N_q, N_b, D]\) where:

  • \(N_f\) = number of facets

  • \(N_q\) = number of quadrature points per facet

  • \(N_b\) = number of basis functions

  • \(D\) = spatial dimension

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Triangular facet gradients of shape \([N_{tf}, N_{tq}, N_b, D]\)

  • Quadrilateral facet gradients of shape \([N_{qf}, N_{qq}, N_b, D]\)

  where:

  • \(N_{tf}\) = number of triangular facets

  • \(N_{tq}\) = quadrature points per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets

  • \(N_{qq}\) = quadrature points per quadrilateral facet

  • \(N_b\) = number of basis functions

  • \(D\) = spatial dimension

property facet_jacobian: Tensor | Tuple[Tensor, Tensor]

Get Jacobian matrices at facet quadrature points.

The facet Jacobian matrix maps derivatives from reference facet coordinates to physical coordinates. For mixed facet elements, returns separate Jacobians for triangular and quadrilateral facets.

The Jacobian matrix \(\mathbf{J}_f = \frac{\partial\mathbf{x}}{\partial\boldsymbol{\xi}_f}\) where \(\mathbf{x}\) are physical coordinates and \(\boldsymbol{\xi}_f\) are reference facet coordinates.

Examples

# For tetrahedral mesh (triangular facets)
transform = Transformation(points, elements, "tetra")
J = transform.facet_jacobian
# J[e,f,q] is Jacobian at qth quadrature point
# on fth facet of eth element

# For prism mesh (mixed facets)
transform = Transformation(points, elements, "wedge")
J_tri, J_quad = transform.facet_jacobian

Returns:

• If element has uniform facets –

  torch.Tensor

  Facet Jacobian tensor of shape \([N_e, N_f, N_q, D-1, D]\) where:

  • \(N_e\) = number of elements

  • \(N_f\) = number of facets per element

  • \(N_q\) = number of quadrature points per facet

  • \(D\) = spatial dimension

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Triangular facet Jacobians of shape \([N_e, N_{tf}, N_{tq}, D-1, D]\)

  • Quadrilateral facet Jacobians of shape \([N_e, N_{qf}, N_{qq}, D-1, D]\)

  where:

  • \(N_e\) = number of elements

  • \(N_{tf}\) = number of triangular facets per element

  • \(N_{tq}\) = quadrature points per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets per element

  • \(N_{qq}\) = quadrature points per quadrilateral facet

  • \(D\) = spatial dimension

property nanson_scale: Tensor | Tuple[Tensor, Tensor]

Get Nanson scale factors for element facets.

The Nanson scale factor \(\mathbf{n}\) is used to transform area elements from reference to physical space:

\[\mathbf{n} = \det(\mathbf{J}) \mathbf{J}^{-T} \mathbf{N} w\]

where:

• \(\mathbf{J}\) is the facet-to-cell Jacobian matrix

• \(\mathbf{N}\) is the outward unit normal in reference space

• \(w\) is the quadrature weight

Examples

# For elements with uniform facets (e.g. tetrahedra)
transform = Transformation(points, elements, "tetra")
n = transform.nanson_scale  # Single tensor for all facets

# For elements with mixed facets (e.g. prisms)
transform = Transformation(points, elements, "wedge")
n_tri, n_quad = transform.nanson_scale  # Separate tensors for tri/quad facets

Returns:

• If element has uniform facets –

  torch.Tensor

  Nanson scale tensor of shape \([N_e, N_f, N_q]\) where:

  • \(N_e\) = number of elements

  • \(N_f\) = number of facets per element

  • \(N_q\) = number of quadrature points per facet

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Triangular facet scales of shape \([N_e, N_{tf}, N_{tq}]\)

  • Quadrilateral facet scales of shape \([N_e, N_{qf}, N_{qq}]\)

  where:

  • \(N_e\) = number of elements

  • \(N_{tf}\) = number of triangular facets per element

  • \(N_{tq}\) = quadrature points per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets per element

  • \(N_{qq}\) = quadrature points per quadrilateral facet

References

Wikiversity, Continuum mechanics — Volume change and area change: https://en.wikiversity.org/wiki/Continuum_mechanics/Volume_change_and_area_change

property phi: Tensor

Get shape function values at quadrature points.

The shape functions \(N_i(\boldsymbol{\xi})\) are evaluated at the quadrature points in the reference element.

Examples

transform = Transformation(points, elements, "triangle")
N = transform.shape_val
# N[i,j] is value of jth shape function at ith quadrature point

Returns:

Shape function values tensor of shape \([N_q, N_b]\) where:

• \(N_q\) = number of quadrature points

• \(N_b\) = number of basis functions

Return type:

Tensor

property gradphi: Tensor

Get shape function gradients at quadrature points.

The shape function gradients \(\nabla N_i(\mathbf{x})\) are evaluated at the quadrature points in physical coordinates.

Examples

transform = Transformation(points, elements, "triangle")
dN = transform.shape_grad
# dN[e,q,i,d] is derivative of ith shape function w.r.t.
# dth coordinate at qth quadrature point in eth element

Returns:

Shape function gradients tensor of shape \([N_e, N_q, N_b, D]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

• \(N_b\) = number of basis functions

• \(D\) = spatial dimension

Return type:

Tensor

property J: Tensor

Get Jacobian matrices at basis points,

The Jacobian matrix \(\mathbf{J} = \frac{\partial\mathbf{x}}{\partial\boldsymbol{\xi}}\) maps derivatives from reference to physical coordinates.

The Jacobian matrix maps derivatives between reference and physical coordinates:

\[\frac{\partial}{\partial x_i} = \sum_{j=1}^D J^{-1}_{ij} \frac{\partial}{\partial \xi_j}\]

where \(J_{ij} = \frac{\partial x_i}{\partial \xi_j}\) and \(D\) is the spatial dimension.

Examples

transform = Transformation(points, elements, "triangle")
J = transform.jacobian
# J[e,q,:,:] is Jacobian matrix at qth quadrature point
# in eth element

Returns:

Jacobian matrices tensor of shape \([N_e, N_q, D, D]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

• \(D\) = spatial dimension

Return type:

Tensor

property detJ: Tensor

Get the determinant of the Jacobian matrix at quadrature points.

The determinant \(\lvert J \rvert = \det\!\left(\frac{\partial\mathbf{x}}{\partial\boldsymbol{\xi}}\right)\) is the local scaling factor between reference and physical coordinates.

Examples

transform = Transformation(points, elements, "triangle")
detJ = transform.detJ
# detJ[e,q] is Jacobian determinant at qth quadrature point
# in eth element

Returns:

Jacobian determinant tensor of shape \([N_e, N_q]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

Return type:

Tensor

property JxW: Tensor

Get the Jacobian determinant times quadrature weights.

The JxW term \(\vert J\vert w\) combines the Jacobian determinant with quadrature weights for numerical integration over physical elements:

\[\int_\Omega f(\mathbf{x}) d\mathbf{x} = \sum_{e=1}^{N_e} \sum_{q=1}^{N_q} f(\mathbf{x}_q^e) \vert J_q^e \vert w_q\]

Examples

transform = Transformation(points, elements, "triangle")
jxw = transform.JxW
# Integrate function over mesh:
integral = torch.sum(f * jxw)  # f evaluated at quadrature points

Returns:

JxW tensor of shape \([N_e, N_q]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

Return type:

Tensor

property G: Tensor

Get Jacobian matrices at basis points,

The Jacobian matrix \(\mathbf{J} = \frac{\partial\mathbf{x}}{\partial\boldsymbol{\xi}}\) maps derivatives from reference to physical coordinates.

The Jacobian matrix maps derivatives between reference and physical coordinates:

\[\frac{\partial}{\partial x_i} = \sum_{j=1}^D J^{-1}_{ij} \frac{\partial}{\partial \xi_j}\]

where \(J_{ij} = \frac{\partial x_i}{\partial \xi_j}\) and \(D\) is the spatial dimension.

Examples

transform = Transformation(points, elements, "triangle")
J = transform.jacobian
# J[e,q,:,:] is Jacobian matrix at qth quadrature point
# in eth element

Returns:

Jacobian matrices tensor of shape \([N_e, N_q, D, D]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

• \(D\) = spatial dimension

Return type:

Tensor

property detG: Tensor

Get the determinant of the Jacobian matrix at quadrature points.

The determinant \(\lvert J \rvert = \det\!\left(\frac{\partial\mathbf{x}}{\partial\boldsymbol{\xi}}\right)\) is the local scaling factor between reference and physical coordinates.

Examples

transform = Transformation(points, elements, "triangle")
detJ = transform.detJ
# detJ[e,q] is Jacobian determinant at qth quadrature point
# in eth element

Returns:

Jacobian determinant tensor of shape \([N_e, N_q]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

Return type:

Tensor

property GxW: Tensor

Get the Jacobian determinant times quadrature weights.

The JxW term \(\vert J\vert w\) combines the Jacobian determinant with quadrature weights for numerical integration over physical elements:

\[\int_\Omega f(\mathbf{x}) d\mathbf{x} = \sum_{e=1}^{N_e} \sum_{q=1}^{N_q} f(\mathbf{x}_q^e) \vert J_q^e \vert w_q\]

Examples

transform = Transformation(points, elements, "triangle")
jxw = transform.JxW
# Integrate function over mesh:
integral = torch.sum(f * jxw)  # f evaluated at quadrature points

Returns:

JxW tensor of shape \([N_e, N_q]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

Return type:

Tensor

property F: Tensor | Tuple[Tensor, Tensor]

Get Jacobian matrices at facet quadrature points.

The facet Jacobian matrix maps derivatives from reference facet coordinates to physical coordinates. For mixed facet elements, returns separate Jacobians for triangular and quadrilateral facets.

The Jacobian matrix \(\mathbf{J}_f = \frac{\partial\mathbf{x}}{\partial\boldsymbol{\xi}_f}\) where \(\mathbf{x}\) are physical coordinates and \(\boldsymbol{\xi}_f\) are reference facet coordinates.

Examples

# For tetrahedral mesh (triangular facets)
transform = Transformation(points, elements, "tetra")
J = transform.facet_jacobian
# J[e,f,q] is Jacobian at qth quadrature point
# on fth facet of eth element

# For prism mesh (mixed facets)
transform = Transformation(points, elements, "wedge")
J_tri, J_quad = transform.facet_jacobian

Returns:

• If element has uniform facets –

  torch.Tensor

  Facet Jacobian tensor of shape \([N_e, N_f, N_q, D-1, D]\) where:

  • \(N_e\) = number of elements

  • \(N_f\) = number of facets per element

  • \(N_q\) = number of quadrature points per facet

  • \(D\) = spatial dimension

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Triangular facet Jacobians of shape \([N_e, N_{tf}, N_{tq}, D-1, D]\)

  • Quadrilateral facet Jacobians of shape \([N_e, N_{qf}, N_{qq}, D-1, D]\)

  where:

  • \(N_e\) = number of elements

  • \(N_{tf}\) = number of triangular facets per element

  • \(N_{tq}\) = quadrature points per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets per element

  • \(N_{qq}\) = quadrature points per quadrilateral facet

  • \(D\) = spatial dimension

property detF: Tensor | Tuple[Tensor, Tensor]

Get determinant of facet Jacobian matrices.

For a facet with normal \(\mathbf{n}\), the facet Jacobian determinant is:

\[|F| = \sqrt{\det(F^T F)}\]

where \(F\) is the facet Jacobian matrix mapping from reference to physical facets.

Examples

# For tetrahedral mesh (triangular facets)
transform = Transformation(points, elements, "tetra")
detF = transform.detF  # Single tensor for tri facets

# For prism mesh (mixed tri/quad facets)
transform = Transformation(points, elements, "wedge")
tri_detF, quad_detF = transform.detF  # Separate tensors

Returns:

• If element has uniform facets –

  torch.Tensor

  Facet Jacobian determinant tensor of shape \([N_e, N_f, N_q]\) where:

  • \(N_e\) = number of elements

  • \(N_f\) = number of facets per element

  • \(N_q\) = number of quadrature points per facet

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Triangular facet determinants of shape \([N_e, N_{tf}, N_{tq}]\)

  • Quadrilateral facet determinants of shape \([N_e, N_{qf}, N_{qq}]\)

  where:

  • \(N_e\) = number of elements

  • \(N_{tf}\) = number of triangular facets

  • \(N_{tq}\) = number of quadrature points per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets

  • \(N_{qq}\) = number of quadrature points per quadrilateral facet

property FxW: Tensor | Tuple[Tensor, Tensor]

Get facet Jacobian determinant times quadrature weights.

The FxW term \(|F|w\) combines the facet Jacobian determinant with quadrature weights for numerical integration over physical facets:

\[\int_{\Gamma} f(\mathbf{x}) d\Gamma = \sum_{e=1}^{N_e} \sum_{f=1}^{N_f} \sum_{q=1}^{N_q} f(\mathbf{x}_q^{ef}) |F_q^{ef}| w_q\]

Examples

# For tetrahedral mesh (triangular facets)
transform = Transformation(points, elements, "tetra")
fxw = transform.FxW  # Single tensor for tri facets

# For prism mesh (mixed tri/quad facets)
transform = Transformation(points, elements, "wedge")
tri_fxw, quad_fxw = transform.FxW  # Separate tensors

# Integrate function over facets:
integral = torch.sum(f * fxw)  # f evaluated at facet quadrature points

Returns:

• If element has uniform facets –

  torch.Tensor

  FxW tensor of shape \([N_e, N_f, N_q]\) where:

  • \(N_e\) = number of elements

  • \(N_f\) = number of facets per element

  • \(N_q\) = number of quadrature points per facet

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Triangular facet FxW of shape \([N_e, N_{tf}, N_{tq}]\)

  • Quadrilateral facet FxW of shape \([N_e, N_{qf}, N_{qq}]\)

  where:

  • \(N_e\) = number of elements

  • \(N_{tf}\) = number of triangular facets

  • \(N_{tq}\) = number of quadrature points per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets

  • \(N_{qq}\) = number of quadrature points per quadrilateral facet

property n: Tensor | Tuple[Tensor, Tensor]

Get Nanson scale factors for element facets.

The Nanson scale factor \(\mathbf{n}\) is used to transform area elements from reference to physical space:

\[\mathbf{n} = \det(\mathbf{J}) \mathbf{J}^{-T} \mathbf{N} w\]

where:

• \(\mathbf{J}\) is the facet-to-cell Jacobian matrix

• \(\mathbf{N}\) is the outward unit normal in reference space

• \(w\) is the quadrature weight

Examples

# For elements with uniform facets (e.g. tetrahedra)
transform = Transformation(points, elements, "tetra")
n = transform.nanson_scale  # Single tensor for all facets

# For elements with mixed facets (e.g. prisms)
transform = Transformation(points, elements, "wedge")
n_tri, n_quad = transform.nanson_scale  # Separate tensors for tri/quad facets

Returns:

• If element has uniform facets –

  torch.Tensor

  Nanson scale tensor of shape \([N_e, N_f, N_q]\) where:

  • \(N_e\) = number of elements

  • \(N_f\) = number of facets per element

  • \(N_q\) = number of quadrature points per facet

• If element has mixed facets –

  Tuple[torch.Tensor, torch.Tensor]

  • Triangular facet scales of shape \([N_e, N_{tf}, N_{tq}]\)

  • Quadrilateral facet scales of shape \([N_e, N_{qf}, N_{qq}]\)

  where:

  • \(N_e\) = number of elements

  • \(N_{tf}\) = number of triangular facets per element

  • \(N_{tq}\) = quadrature points per triangular facet

  • \(N_{qf}\) = number of quadrilateral facets per element

  • \(N_{qq}\) = quadrature points per quadrilateral facet

References

Wikiversity, Continuum mechanics — Volume change and area change: https://en.wikiversity.org/wiki/Continuum_mechanics/Volume_change_and_area_change

batch_quadrature(start: int, batch: int) → Tuple[Tensor, Tensor]

Get a batch of quadrature points and weights.

Examples

transform = Transformation(points, elements, "triangle")

# Get first 10 quadrature points
w, q = transform.batch_quadrature(start=0, batch=10)

# Get all quadrature points after index 5
w, q = transform.batch_quadrature(start=5, batch=-1)

Parameters:

• start (int) – The starting index of the batch

• batch (int) – The number of quadrature points in the batch. If -1, returns all points starting from start

Returns:

• Quadrature weights tensor of shape \([N_b]\)

• Quadrature points tensor of shape \([N_b, D]\)

where:

• \(N_b\) = batch size

• \(D\) = spatial dimension

Return type:

Tuple[Tensor, Tensor]

batch_elements_coords(start: int, batch: int) → Tensor

Get physical coordinates for a batch of elements.

Examples

transform = Transformation(points, elements, "triangle")

# Get coordinates for first 5 elements
coords = transform.batch_elements_coords(start=0, batch=5)

# Get coordinates for all elements after index 10
coords = transform.batch_elements_coords(start=10, batch=-1)

Parameters:

• start (int) – The starting element index

• batch (int) – The number of elements in the batch. If -1, returns all elements starting from start

Returns:

Element coordinates tensor of shape \([N_b, N_n, D]\) where:

• \(N_b\) = batch size

• \(N_n\) = nodes per element

• \(D\) = spatial dimension

Return type:

Tensor

batch_shape_val(start: int, batch: int) → Tensor

Get shape function values for a batch of quadrature points.

Examples

transform = Transformation(points, elements, "triangle")

# Get shape values at first 10 quadrature points
vals = transform.batch_shape_val(start=0, batch=10)

# Get shape values at all points after index 5
vals = transform.batch_shape_val(start=5, batch=-1)

Parameters:

• start (int) – The starting quadrature point index

• batch (int) – The number of quadrature points in the batch

Returns:

Shape function values tensor of shape \([N_b, N_s]\) where:

• \(N_b\) = batch size

• \(N_s\) = shape functions per element

Return type:

Tensor

batch_shape_grad_jxw(element_start: int = 0, element_batch: int = -1, quadrature_start: int = 0, quadrature_batch: int = -1) → Tuple[Tensor, Tensor]

Get shape function gradients and Jacobian weights for batches of elements and quadrature points.

Examples

transform = Transformation(points, elements, "triangle")

# Get values for first 5 elements and first 10 quadrature points
grads, jxw = transform.batch_shape_grad_jxw(
    element_start=0, element_batch=5,
    quadrature_start=0, quadrature_batch=10
)

# Get values for all elements and quadrature points
grads, jxw = transform.batch_shape_grad_jxw()

Parameters:

• element_start (int, optional) – The starting element index, by default 0

• element_batch (int, optional) – The number of elements in the batch. If -1, uses all elements after start, by default -1

• quadrature_start (int, optional) – The starting quadrature point index, by default 0

• quadrature_batch (int, optional) – The number of quadrature points in the batch. If -1, uses all points after start, by default -1

Returns:

• Shape gradients tensor of shape \([N_e, N_q, N_s, D]\)

• Jacobian weights tensor of shape \([N_e, N_q]\)

where:

• \(N_e\) = number of elements in batch

• \(N_q\) = number of quadrature points in batch

• \(N_s\) = shape functions per element

• \(D\) = spatial dimension

Return type:

Tuple[Tensor, Tensor]

Reference-shape catalog

TensorMesh ships seven reference shapes. The class objects double as types (you name them in reorder() calls, in element_type2element() results, as the value of element, …) and as a small public surface for shape introspection — basis-node coordinates, quadrature rules, facet type, Gmsh/VTK permutation. Users rarely call methods on subclasses directly; the methods listed below are the ones referenced from the user guide and the example gallery, plus the extension hooks an element-type author overrides.

class Element

Bases: object

Base class for the seven reference shapes (line, triangle, quad, tet, hex, pyramid, prism).

Each subclass is a thin namespace of class methods — get_basis, get_quadrature, get_facet_type, reorder, … — that the assembler and mesh I/O paths call internally to obtain interpolation-node layouts, shape functions, quadrature rules, and the Gmsh/VTK ↔ TensorMesh permutation.

Subclasses are not instantiated; Triangle, Hexahedron, etc. are used directly as types. See Elements and Quadrature for the full API tour and Gmsh / VTK ↔ TensorMesh node ordering for a visual comparison of the node-numbering conventions.

points: Tensor

Coordinates of element vertices Shape: \([V, D]\) where \(V\) = number of vertices and \(D\) = spatial dimension

vertex: Tensor

Vertex indices for the element Shape: \([V, 1]\) where \(V\) = number of vertices

edge: Tensor

Edge connectivity defining pairs of vertex indices Shape: \([E, 2]\) where \(E\) = number of edges

face: Tuple[Tuple[int, ...], ...] | None

Face connectivity defining vertex indices for each face. Shape: Tuple of tuples, each inner tuple contains vertex indices for one face.

Example:

((0,1,2), (1,2,3)) for two triangular faces

cell: Tensor | None

Cell connectivity defining vertex indices for volumetric elements Shape: \([N_c, V_c]\) where \(N_c\) = number of cells and \(V_c\) = vertices per cell

dim: int

Spatial dimension of the element (1D, 2D, or 3D)

n_vertex: int

Number of vertices in the element

n_edge: int

Number of edges in the element

n_face: int

Number of faces in the element

n_cell: int

Number of cells in the element

is_mix_facet: bool = False

Boolean indicating whether the element has mixed facet types.

False indicates uniform facets (e.g. all triangles or all quads). True indicates mixed facets (e.g. prism or pyramid).

classmethod reorder(elements: Tensor, to_gmsh: bool = True) → Tensor

Reorder element connectivity between TensorMesh internal ordering and Gmsh/VTK ordering.

TensorMesh stores connectivity for some tensor-product elements (e.g. Quad, Hex) in a lexicographic node order that matches the basis / tensor-product polynomial layout. However, file formats like VTK/VTU (and Gmsh conventions) expect a different ordering.

This helper provides a single, centralized conversion point so that: - internal FEM kernels always see TensorMesh ordering - IO / visualization (VTK/VTU) can export in Gmsh/VTK ordering

Parameters:

• elements – Connectivity tensor of shape [..., n_nodes] (typically [n_elem, n_nodes]).

• to_gmsh –

  • False: convert Gmsh/VTK -> TensorMesh (used when reading meshes)

  • True: convert TensorMesh -> Gmsh/VTK (used when writing VTK/VTU)

Returns:

Reordered connectivity tensor with the same shape/dtype/device as input.

Return type:

Tensor

classmethod get_gmsh_permutation(n_nodes: int, device: device = torch.device('cpu')) → Tensor

Return a permutation perm (shape [n_nodes]) mapping Gmsh/VTK -> TensorMesh.

Given an input connectivity conn_gmsh[..., :] in Gmsh/VTK ordering, the internal ordering is:

conn_internal = conn_gmsh.index_select(-1, perm)

Notes

Subclasses should override this for element types where TensorMesh ordering differs from the external Gmsh/VTK ordering. If not overridden, the default is identity.

classmethod get_facet_type() → Type[Element] | Tuple[Type[Element], Type[Element]]

Get the element type(s) for the facets of this element.

This method returns the element type(s) that make up the facets (faces) of the current element. For example, a triangle element has line segment facets, while a prism has both triangular and quadrilateral facets.

Examples

>>> from tensormesh.element import Triangle, Line
>>> Triangle.get_facet_type()
Line

>>> from tensormesh.element import Prism, Triangle, Quadrilateral
>>> facet_types = Prism.get_facet_type()
>>> facet_types == (Triangle, Quadrilateral)
True

>>> from tensormesh.element import Tetrahedron
>>> Tetrahedron.get_facet_type()
Triangle

Returns:

The element type(s) for the facets. Returns either:

• A single Element type for elements with uniform facets (e.g. Line for Triangle)

• A tuple of Element types for elements with mixed facets (e.g. (Triangle, Quadrilateral) for Prism)

Return type:

Union[Type[’Element’], Tuple[Type[’Element’], Type[’Element’]]]

classmethod get_basis(order: int = 1, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Tensor

Get the basis node coordinates for the element.

This method returns the coordinates of the basis nodes that define the element’s shape functions. The number and location of basis nodes depends on the polynomial order.

For example: - For order=1 (linear), returns vertices of the element - For order=2 (quadratic), returns vertices plus midpoints of edges - For higher orders, returns additional nodes in the interior

The basis nodes are used to construct the element’s shape functions and define the geometric mapping between reference and physical elements.

Examples

>>> from tensormesh.element import Triangle
>>> basis = Triangle.get_basis(order=1)
>>> basis.shape
torch.Size([3, 2])  # 3 vertices in 2D
>>> basis
tensor([[0., 0.],
        [1., 0.],
        [0., 1.]])

>>> # Quadratic basis has additional nodes at edge midpoints
>>> basis = Triangle.get_basis(order=2)
>>> basis.shape
torch.Size([6, 2])  # 3 vertices + 3 edge midpoints

>>> from tensormesh.element import Tetrahedron
>>> basis = Tetrahedron.get_basis(order=1)
>>> basis.shape
torch.Size([4, 3])  # 4 vertices in 3D
>>> basis
tensor([[0., 0., 0.],
        [1., 0., 0.],
        [0., 1., 0.],
        [0., 0., 1.]])

Parameters:

• order (int) – the order of the basis

• dtype (dtype) – the float data type of the polynomial

• device (device) – the device of the polynomial

Returns:

basis – 2D Tensor of shape [B, D], where B is the number of basis nodes and D is the dimension

Return type:

Tensor

classmethod get_facet(order: int = 1) → Tensor | Tuple[Tensor, Tensor]

Get the facet connectivity for the element.

For elements with uniform facets (e.g. Triangle, Quadrilateral), returns a single tensor containing facet connectivity. For elements with mixed facets (e.g. Prism), returns a tuple of tensors containing facet connectivity for each facet type.

The facet connectivity defines how the basis nodes are connected to form each facet of the element. For example, for a triangle element with linear basis (order=1), the facet connectivity would be:

[[0,1],  # First edge connects nodes 0 and 1
 [1,2],  # Second edge connects nodes 1 and 2
 [2,0]]  # Third edge connects nodes 2 and 0

For a prism element which has both triangular and quadrilateral facets, would return:

(tri_facets, quad_facets)

Where tri_facets contains connectivity for the triangular faces and quad_facets contains connectivity for the quadrilateral faces.

Examples

>>> from tensormesh.element import Triangle
>>> facets = Triangle.get_facet(order=1)
>>> facets
tensor([[0, 1],  # First edge connects nodes 0 and 1
        [1, 2],  # Second edge connects nodes 1 and 2
        [2, 0]]) # Third edge connects nodes 2 and 0

>>> # Higher order elements have more nodes per facet
>>> facets = Triangle.get_facet(order=2)
>>> facets.shape
torch.Size([3, 3])  # 3 edges with 3 nodes each

>>> from tensormesh.element import Prism
>>> tri_facets, quad_facets = Prism.get_facet(order=1)
>>> tri_facets.shape  # Two triangular faces
torch.Size([2, 3])
>>> quad_facets.shape  # Three quadrilateral faces
torch.Size([3, 4])

Parameters:

order (int) – the order of the basis

Returns:

For elements with uniform facets, a 2D tensor of shape \([F, B_f]\) where \(F\) is the number of facets and \(B_f\) the number of basis nodes per facet.

For elements with mixed facets, a pair (tri_facet, quad_facet) of shapes \([F_t, B_{ft}]\) and \([F_q, B_{fq}]\), respectively.

Return type:

Tensor or Tuple[Tensor, Tensor]

classmethod get_edge(order: int) → Tensor

Get the edge basis indices for the element.

Returns the connectivity array for the edges of the element. Each row represents an edge and contains the indices of the basis nodes along that edge.

Examples

>>> from tensormesh.element import Triangle
>>> edges = Triangle.get_edge(order=1)
>>> edges
tensor([[0, 1],  # First edge connects nodes 0 and 1
        [1, 2],  # Second edge connects nodes 1 and 2
        [2, 0]]) # Third edge connects nodes 2 and 0

>>> # Higher order elements have more nodes per edge
>>> edges = Triangle.get_edge(order=2)
>>> edges.shape
torch.Size([3, 3])  # 3 edges with 3 nodes each

Parameters:

order (int) – the order of the basis

Returns:

edge – 2D Tensor of shape [n_edge, 2] containing edge connectivity

Return type:

Tensor

classmethod get_polynomial(order: int = 1, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomial

Get the polynomial form for the element.

For a given order n, returns a polynomial with terms up to maximum total degree n. The polynomial will have n_vars equal to the dimension of the element.

For example,

• For a 2D element with order=1: \(P(x,y) = a_0 + a_1x + a_2y\)

• For order=2: \(P(x,y) = a_0 + a_1x + a_2y + a_3x^2 + a_4xy + a_5y^2\)

The exact form and number of terms depends on the specific element type.

Examples

>>> from tensormesh.element import Triangle
>>> poly = Triangle.get_polynomial(order=1)
>>> poly.n_vars  # 2D element
2
>>> poly.n_terms  # Linear terms: 1, x, y
3

>>> poly = Triangle.get_polynomial(order=2)
>>> poly.n_terms  # Quadratic terms: 1, x, y, x^2, xy, y^2
6

>>> from tensormesh.element import Tetrahedron
>>> poly = Tetrahedron.get_polynomial(order=1)
>>> poly.n_vars  # 3D element
3
>>> poly.n_terms  # Linear terms: 1, x, y, z
4

Parameters:

• order (int) – the order of the basis

• dtype (dtype) – the float data type of the polynomial

• device (device) – the device of the polynomial

Returns:

polynomial – A polynomial object with number of variables equal to the element dimension,

Return type:

Polynomial

classmethod get_quadrature(order: int = 1, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Tuple[Tensor, Tensor]

Get the quadrature points and weights for the element.

For elements with single type of facets (e.g. Tetrahedron):

• Returns a tuple of two tensors containing quadrature weights and points

For elements with mixed facets (e.g. Prism):

• Returns a tuple of four tensors containing quadrature weights and points for each facet type

Examples

>>> from tensormesh.element import Triangle
>>> weights, points = Triangle.get_quadrature(order=1)
>>> weights.shape  # Number of quadrature points
torch.Size([1])
>>> points.shape  # [n_points, dimension]
torch.Size([1, 2])

>>> weights, points = Triangle.get_quadrature(order=2)
>>> weights.shape  # More quadrature points for higher order
torch.Size([3])
>>> points.shape
torch.Size([3, 2])

>>> from tensormesh.element import Tetrahedron
>>> weights, points = Tetrahedron.get_quadrature(order=1)
>>> weights.shape
torch.Size([1])
>>> points.shape  # 3D points
torch.Size([1, 3])

Parameters:

• order (int) – quadrature order

• dtype (dtype) – the float data type of the polynomial

• device (device) – the device of the polynomial

Returns:

• quadrature_weights (torch.Tensor) – 1D tensor of shape \([N_q]\) where \(N_q\) is the number of quadrature points.

• quadrature_points (torch.Tensor) – 2D tensor of shape \([N_q, D]\) where \(D\) is the element dimension (1, 2, or 3).

classmethod get_facet_quadrature(order: int = 1, transform: bool = True, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Tuple[Tensor, Tensor] | Tuple[Tensor, Tensor, Tensor, Tensor]

Get the quadrature points and weights for the facets of the element.

This method returns quadrature rules for integrating over element facets. The return format depends on whether the element has uniform or mixed facet types.

For elements with uniform facets (e.g. Tetrahedron):

• Returns a tuple of (weights, points) tensors

• weights: quadrature weights for each facet point

• points: quadrature point coordinates

For elements with mixed facets (e.g. Prism):

• Returns a tuple of (tri_weights, tri_points, quad_weights, quad_points)

• tri_weights: quadrature weights for triangular facets

• tri_points: quadrature points for triangular facets

• quad_weights: quadrature weights for quadrilateral facets

• quad_points: quadrature points for quadrilateral facets

The quadrature points can be returned either in the reference element coordinates (transform=True) or in facet-local coordinates (transform=False).

Examples

>>> # Example for tetrahedron element with uniform triangular facets
>>> weights, points = Tetrahedron.get_facet_quadrature(order=1)
>>> weights.shape  # [n_facet=4, n_quad_points=1]
torch.Size([4, 1])
>>> points.shape  # [n_facet=4, n_quad_points=1, dim=3]
torch.Size([4, 1, 3])

>>> # Example for prism element with mixed facets (triangular and quadrilateral)
>>> tri_w, tri_p, quad_w, quad_p = Prism.get_facet_quadrature(order=1)
>>> tri_w.shape  # [n_tri_facet=2, n_quad_points=1]
torch.Size([2, 1])
>>> tri_p.shape  # [n_tri_facet=2, n_quad_points=1, dim=3]
torch.Size([2, 1, 3])
>>> quad_w.shape # [n_quad_facet=3, n_quad_points=1]
torch.Size([3, 1])
>>> quad_p.shape # [n_quad_facet=3, n_quad_points=1, dim=3]
torch.Size([3, 1, 3])

>>> # Example with transform=False returns points in facet-local coordinates
>>> weights, points = Triangle.get_facet_quadrature(order=1, transform=False)
>>> weights.shape  # [n_quad_points=1]
torch.Size([1])
>>> points.shape  # [n_quad_points=1, dim-1=1]
torch.Size([1, 1])

Parameters:

• order (int, optional) – Quadrature order. Default is 1.

• transform (bool, optional) – Whether to return the transformed facet quadrature. If True, returns points in reference element coordinates with shape \([N_f, N_q, D]\), where \(N_f\) = number of facets, \(N_q\) = quadrature points per facet, and \(D\) = spatial dimension. If False, returns points in facet-local coordinates with shape \([N_q, D-1]\), where \(N_q\) = quadrature points per facet and \(D\) = spatial dimension.

• dtype (dtype, optional) – The float data type of the polynomial.

• device (device, optional) – The device of the polynomial.

Returns:

For elements with uniform facets, a pair (weights, points) whose shapes depend on transform:

• transform=True — weights of shape \([N_f, N_q]\), points of shape \([N_f, N_q, D]\) (cell coordinates).

• transform=False — weights of shape \([N_q]\), points of shape \([N_q, D-1]\) (facet-local coordinates).

For elements with mixed facets, a 4-tuple (tri_weights, tri_points, quad_weights, quad_points) with the obvious analogues of the shapes above for the triangular and quadrilateral facet kinds.

Return type:

Tuple[Tensor, Tensor] or Tuple[Tensor, Tensor, Tensor, Tensor]

classmethod get_basis_fns(order: int = 1, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomials

Get the basis functions for this element type.

Returns polynomial basis functions that interpolate values at the element nodes. The basis functions have the Kronecker delta property, meaning each basis function is 1 at its corresponding node and 0 at all other nodes.

For example, for a linear triangle element (order=1), returns 3 basis functions corresponding to the 3 vertices. Each basis function is 1 at one vertex and 0 at the other vertices.

For higher orders, returns basis functions for both vertex and edge nodes. The number and location of nodes depends on the specific element type.

The basis functions are returned as a Polynomials object containing the polynomial coefficients and exponents. This allows efficient evaluation and manipulation of the basis functions.

Examples

>>> from tensormesh.element import Triangle
>>> basis_fns = Triangle.get_basis_fns(order=1)
>>> basis_fns
Polynomials([
    1 - x - y,  # First basis function, 1 at vertex 0 (0,0)
    x,          # Second basis function, 1 at vertex 1 (1,0)
    y           # Third basis function, 1 at vertex 2 (0,1)
])

>>> # Evaluate basis functions at a point
>>> point = torch.tensor([[0.5, 0.25]])
>>> values = basis_fns(point)
>>> values
tensor([[0.25, 0.50, 0.25]])  # Values sum to 1

>>> # Higher order basis functions
>>> basis_fns = Triangle.get_basis_fns(order=2)
>>> basis_fns
Polynomials([
    # Vertex basis functions
    2*x^2 + 2*y^2 - 3*x - 3*y + 1,  # Node 0
    2*x^2 - x,                       # Node 1
    2*y^2 - y,                       # Node 2
    # Edge basis functions
    4*x*y,                           # Node 3
    -4*x*y + 4*x,                    # Node 4
    -4*x*y + 4*y                     # Node 5
])

>>> from tensormesh.element import Tetrahedron
>>> basis_fns = Tetrahedron.get_basis_fns(order=1)
>>> basis_fns
Polynomials([
    1 - x - y - z,  # First basis function
    x,              # Second basis function
    y,              # Third basis function
    z               # Fourth basis function
])

Parameters:

• order (int) – The order of the basis functions.

• dtype (dtype) – The float data type of the polynomial.

• device (device) – The device to place the polynomial on.

Returns:

The basis functions of the element. Shape is \([N_b]\) where \(N_b\) is the number of basis functions, with \(N_{vars}=D\) variables and \(N_{exp}=N_b\) expansion terms.

Return type:

Polynomials

classmethod get_basis_grad_fns(order: int = 1, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomials

Get the gradient of basis functions with respect to reference coordinates.

For a basis function \(\phi(x_1,\ldots,x_n)\), computes the gradient:

\[\nabla \phi = \left[\frac{\partial \phi}{\partial x_1}, \ldots, \frac{\partial \phi}{\partial x_n}\right]\]

The gradient is computed by taking the partial derivative of each basis function with respect to each coordinate direction \(x_i\) in the reference element.

Examples

>>> # Get gradients of linear basis functions for a triangle
>>> from tensormesh.element import Triangle
>>> grad_fns = Triangle.get_basis_grad_fns(order=1)
>>> points = torch.tensor([[0.0, 0.0], [0.5, 0.5]])
>>> grads = grad_fns(points)  # [n_points=2, dim=2, n_basis=3]
>>> grads.shape
torch.Size([2, 2, 3])

>>> # Get gradients of quadratic basis functions for a tetrahedron
>>> from tensormesh.element import Tetrahedron
>>> grad_fns = Tetrahedron.get_basis_grad_fns(order=2)
>>> points = torch.tensor([[0.25, 0.25, 0.25]])
>>> grads = grad_fns(points)  # [n_points=1, dim=3, n_basis=10]
>>> grads.shape
torch.Size([1, 3, 10])

Parameters:

• order (int) – The order of the basis functions

• dtype (dtype) – The float data type of the polynomial

• device (device) – The device to place the polynomial on

Returns:

The gradient basis functions. Shape is \([D, N_b]\) where - \(D\) = spatial dimension - \(N_b\) = number of basis functions - n_vars = \(D\) - n_terms = \(N_b\)

Return type:

Polynomials

classmethod eval_cell_jacobian(quadrature: Tensor, element_coords: Tensor, basis_order: int = 1) → Tensor

Evaluate the Jacobian matrix at quadrature points for each element.

The Jacobian matrix represents the linear transformation from reference coordinates to physical coordinates. It is computed by multiplying the element coordinates by the gradients of the basis functions.

For each quadrature point q and element e, computes:

\[J_{e,q,i,j} = \sum_b \text{coords}_{e,b,j} \frac{\partial \phi_b}{\partial x_i}(q)\]

where:

• \(\text{coords}_{e,b,j}\) is the j-th coordinate of basis node b in element e

• \(\phi_b\) is basis function b

• \(\frac{\partial \phi_b}{\partial x_i}\) is the derivative of basis function b with respect to reference coordinate i

Examples

element = Triangle
quadrature = torch.tensor([[0.0, 0.0], [0.5, 0.5]], dtype=torch.float32)
element_coords = torch.tensor([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]], dtype=torch.float32)
jacobian = element.eval_cell_jacobian(quadrature, element_coords)
print(jacobian.shape)  # torch.Size([1, 2, 2, 2])

Parameters:

• quadrature (Tensor) –

  Shape \([N_q, D]\) or \([N_e, N_q, D]\) where:

  • \(N_q\) = number of quadrature points

  • \(N_e\) = number of elements

  • \(D\) = spatial dimension

  Quadrature points for each element

• element_coords (Tensor) –

  Shape \([N_e, N_b, D]\) where:

  • \(N_e\) = number of elements

  • \(N_b\) = number of basis nodes

  • \(D\) = spatial dimension

  Physical coordinates of element nodes

• basis_order (int, optional) – Order of basis functions, by default 1

Returns:

Shape \([N_e, N_q, D, D]\) where:

• \(N_e\) = number of elements

• \(N_q\) = number of quadrature points

• \(D\) = spatial dimension

Jacobian matrix evaluated at each quadrature point for each element

Return type:

Tensor

classmethod eval_shape_val(quadrature_points: Tensor | None = None, order: int = 1, quadrature_order: int = 1, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Tensor

Evaluates the shape functions at quadrature points.

For a given set of quadrature points, evaluates each basis function phi_i at those points. Returns a matrix where entry (q,i) is the value of basis function i at quadrature point q.

The shape functions have two key properties:

1. Partition of unity: The functions sum to 1 at any point

2. Interpolation property: phi_i(x_j) = delta_ij at the nodes

For example, for linear elements (order=1):

\[\begin{split}\phi_1(x) &= 1-x \quad \text{(equals 1 at x=0, 0 at x=1)} \\ \phi_2(x) &= x \quad \text{(equals 0 at x=0, 1 at x=1)}\end{split}\]

For quadratic elements (order=2), includes additional mid-node basis functions.

If quadrature_points are not provided, uses default quadrature rule of specified order.

Examples

import torch
from tensormesh.element import Element

# Create quadrature points
quad_points = torch.tensor([[0.0], [0.5], [1.0]])

# Evaluate shape functions at quadrature points
shape_vals = Element.eval_shape_val(quad_points, order=1)
print(shape_vals.shape)  # [3, 2] for linear element

Parameters:

• quadrature_points (Optional[Tensor]) – Tensor of shape \([N_q, D]\) containing quadrature points, where \(N_q\) = number of quadrature points and \(D\) = spatial dimension

• order (int) – Order of the basis functions

• quadrature_order (int) – The order for quadrature, if quadrature is not None, quadrature_order will be ignored

• dtype (dtype) – The float data type of the shape value, if quadrature_points is not None, the dtype will be ignored. Default is torch.float32

• device (device) – The device of the shape value, if quadrature_points is not None, the device will be ignored. Default is torch.device(‘cpu’)

Returns:

Tensor of shape \([N_q, N_b]\) containing shape function values, where \(N_q\) = number of quadrature points and \(N_b\) = number of basis functions

Return type:

Tensor

classmethod eval_shape_grad(element_coords: Tensor, quadrature: Tensor | None = None, basis_order: int = 1, quadrature_order: int = 1) → Tuple[Tensor, Tensor]

Evaluates the shape function gradients and cell Jacobians at quadrature points.

For each element, computes:

1. Shape function gradients in reference coordinates

2. Cell Jacobian mapping from reference to physical coordinates

3. Transforms gradients to physical coordinates using inverse Jacobian

The shape function gradients are used in weak form assembly:

\[\nabla \phi_i = \sum_j J^{-1}_{jk} \frac{\partial \phi_i}{\partial \xi_k}\]

where:

• \(\phi_i\) are the basis functions

• \(J\) is the Jacobian matrix

• \(\xi_k\) are reference coordinates

The cell Jacobian \(J_{ij} = \frac{\partial x_i}{\partial \xi_j}\) maps derivatives between reference and physical coordinates.

For linear elements (order=1), the Jacobian is constant within each element. For higher order elements, the Jacobian varies spatially.

Examples

import torch
from tensormesh.element import Element

# Create element coordinates
element_coords = torch.tensor([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]])

# Evaluate shape function gradients and Jacobian
shape_grad, cell_jacobian = Element.eval_shape_grad(element_coords, basis_order=1)
print(shape_grad.shape)      # [1, n_quad, 3, 2]
print(cell_jacobian.shape)   # [1, n_quad, 2, 2]

Parameters:

• element_coords (Tensor) – Tensor of shape \([N_e, N_b, D]\) containing element vertex coordinates, where: - \(N_e\) = number of elements - \(N_b\) = number of basis functions - \(D\) = spatial dimension

• quadrature (Tensor) – Tensor of shape \([N_q, D]\) containing quadrature points for each element, where: - \(N_q\) = number of quadrature points - \(D\) = spatial dimension

• basis_order (int) – Order for basis functions

• quadrature_order (int) – Order for quadrature, if element_coords is not None, the quadrature_order will be ignored, default is 1

Returns:

• shape_grad (torch.Tensor) – Tensor of shape \([N_e, N_q, N_b, D]\) containing shape function gradients, where: - \(N_e\) = number of elements - \(N_q\) = number of quadrature points - \(N_b\) = number of basis functions - \(D\) = spatial dimension

• cell_jacobian (torch.Tensor) – Tensor of shape \([N_e, N_q, D, D]\) containing cell Jacobians, where: - \(N_e\) = number of elements - \(N_q\) = number of quadrature points - \(D\) = spatial dimension

classmethod get_local_facet_mapping_fns(dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomials

Maps points from the reference facet to the facet of the reference element.

For a facet \(f\), defines a mapping \(\gamma_f(r)\) from reference facet coordinates \(r\) to reference element coordinates \(\xi\):

\[\gamma_f: \hat{F} \rightarrow F\]

where: - \(\hat{F}\) is the reference facet (e.g. line for 2D, triangle for 3D) - \(F\) is the facet of the reference element - \(\gamma_f(r)\) is a linear mapping

Examples

import torch
from tensormesh.element import Triangle

# Get mapping functions for triangle element facets
facet_fns = Triangle.get_local_facet_mapping_fns()

# Map points from reference facet to element facet
ref_points = torch.tensor([[0.0], [0.5], [1.0]])  # Points on reference line
mapped_points = facet_fns.map(ref_points)  # Points on triangle facet

Parameters:

• dtype (dtype) – the float data type of the polynomials default is torch.float32

• device (device) – the device of the polynomials default is torch.device(‘cpu’)

Returns:

local_facet_fns – Polynomial functions mapping reference facet to element facet coordinates. Shape is \([N_f, D]\) where:

• \(N_f\) = number of facets

• \(D\) = spatial dimension

Each polynomial has:

• \(D-1\) variables (reference facet coordinates)

• \(D\) terms (affine transformation)

The mapping takes the form:

\[\gamma_f(r) = c + \sum_{i=1}^{D-1} \alpha_i r_i\]

where \(c, \alpha_i \in \mathbb{R}^D\) are the coefficients

Return type:

Polynomials

classmethod get_local_facet_mapping_grad_fns(dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomials

Computes the gradient of the facet mapping function with respect to the reference coordinates.

For a facet mapping function \(\gamma_f(r)\) that maps from reference coordinates \(r\) to global coordinates, computes the Jacobian matrix:

\[\frac{\partial \gamma_f (r)}{\partial r}\]

This gives the derivatives of each global coordinate with respect to each reference coordinate.

Examples

element = Triangle()
facet_mapping_grad = element.get_local_facet_mapping_grad_fns()
# Returns polynomial tensor with derivatives of mapping
# Shape: [1, 3, 2] for 2D triangle element (1 ref coord, 3 edges, 2 global coords)

Parameters:

• dtype (dtype) – the float data type of the polynomials default is torch.float32

• device (device) – the device of the polynomials default is torch.device(‘cpu’)

Returns:

Tensor containing the Jacobian with shape \([D-1, N_f, D]\) where:

• \(D-1\) = number of reference coordinates (gradient dimension)

• \(N_f\) = number of facets

• \(D\) = spatial dimension (global coordinates)

The polynomial has \(n_{vars}=D-1\) variables and \(n_{terms}=D\) terms.

Return type:

PolynomialTensor

classmethod get_outwards_facet_normal(dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Tensor

Computes the outward unit normal vectors for each facet of the element.

This returns the geometric normal vectors without the Nanson scale factor. For 2D elements, returns normals for each edge. For 3D elements, returns normals for each face.

Examples

element = Triangle()
normals = element.get_outwards_facet_normal()
# Returns tensor of shape [3, 2] containing outward normals
# for each edge of the triangle

Parameters:

• dtype (dtype) – The float data type for the output tensor Default is torch.float32

• device (device) – The device to place the output tensor on Default is torch.device(‘cpu’)

Returns:

outward_facet_normal – Tensor of shape \([N_f, D]\) containing the outward unit normal vector for each facet, where:

• \(N_f\) = number of facets (edges in 2D, faces in 3D)

• \(D\) = spatial dimension (2 or 3)

Return type:

Tensor

classmethod get_n_facet() → int

Number of facets of the reference element.

Edges for 2D shapes, faces for 3D shapes.

Returns:

Number of facets.

Return type:

int

classmethod get_n_basis(order: int = 1) → int

Number of basis (interpolation) nodes at the given order.

Parameters:

order (int, optional) – Polynomial order. Defaults to 1.

Returns:

Number of basis nodes.

Return type:

int

classmethod eval_facet_cell_jacobian(element_coords: Tensor, basis_order: int = 1, quadrature_order: int = 1) → Tensor | Tuple[Tensor, Tensor]

Evaluates the Jacobian matrix of the mapping from reference to physical coordinates at facet quadrature points.

For each element and facet, computes the Jacobian matrix \(\mathbf{J} = \frac{\partial\mathbf{x}}{\partial\boldsymbol{\xi}}\) at the quadrature points on that facet. The Jacobian describes how the reference element coordinates \(\boldsymbol{\xi}\) map to physical coordinates \(\mathbf{x}\).

For mixed elements with both triangular and quadrilateral facets (e.g. pyramids), returns separate Jacobians for each facet type.

The Jacobian matrix components are:

\[J_{ij} = \frac{\partial x_i}{\partial \xi_j}\]

Examples

# Create a triangle element
element = Triangle()
# Element coordinates for a single triangle
coords = torch.tensor([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]])
# Evaluate facet Jacobians using linear basis and 2nd order quadrature
jacobians = element.eval_facet_cell_jacobian(coords, basis_order=1, quadrature_order=2)
# jacobians shape: [1, 3, 2, 2, 2] for 1 element, 3 edges, 2 quad points per edge

Parameters:

• element_coords (Tensor) – 3D Tensor of shape \([N_e, N_b, D]\) containing the physical coordinates of element nodes, where \(N_e\) = number of elements, \(N_b\) = number of basis nodes, \(D\) = spatial dimension

• basis_order (int) – Order of the basis functions used for the mapping

• quadrature_order (int) – Order of quadrature rule used for facet integration

Returns:

For elements with uniform facets, a 5D tensor of shape \([N_e, N_f, N_q, D, D]\).

For elements with mixed facets (Prism, Pyramid), a pair (tri_facet_cell_jacobian, quad_facet_cell_jacobian) of shapes \([N_e, N_{tf}, N_{tq}, D, D]\) and \([N_e, N_{qf}, N_{qq}, D, D]\).

Return type:

Tensor or Tuple[Tensor, Tensor]

Notes

The Jacobian is used to:

• Transform integrals from reference to physical space:

  \[\int_K f(\mathbf{x}) \,d\mathbf{x} = \int_{\hat{K}} f(\mathbf{x}(\boldsymbol{\xi})) \vert\det(\mathbf{J})\vert \,d\boldsymbol{\xi}\]

• Compute geometric quantities like normal vectors and surface measures

• Map derivatives between reference and physical coordinates:

  \[\nabla_{\mathbf{x}} = \mathbf{J}^{-T} \nabla_{\boldsymbol{\xi}}\]

classmethod eval_facet_jacobian(element_coords: Tensor, basis_order: int = 1, quadrature_order: int = 1) → Tensor | Tuple[Tensor, Tensor]

Evaluate the Jacobian mapping from reference facet coordinates to physical coordinates.

For each facet of the element, computes the Jacobian matrix that maps from local facet coordinates to global physical coordinates.

The facet Jacobian is used to:

• Transform integrals over facets from reference to physical space:

  \[\int_F f(\mathbf{x}) \,ds = \int_{\hat{F}} f(\mathbf{x}(\boldsymbol{\xi})) \|\mathbf{J}\mathbf{t}\| \,d\boldsymbol{\xi}\]

• Compute geometric quantities like normal vectors and surface measures

• Map derivatives between reference facet and physical coordinates:

  \[\nabla_{\mathbf{x}} = \mathbf{J}^{-T} \nabla_{\boldsymbol{\xi}}\]

where \(\mathbf{t}\) is a tangent vector to the facet.

For mixed elements like pyramids that have both triangular and quadrilateral facets, returns separate Jacobians for each facet type.

Examples

# Create a triangle element
element_coords = torch.tensor([[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0]]])
# Evaluate facet Jacobian with linear basis and 2nd order quadrature
facet_jacobian = Triangle.eval_facet_jacobian(element_coords,
                                            basis_order=1,
                                            quadrature_order=2)
facet_jacobian.shape  # torch.Size([1, 3, 2, 1, 2])

Parameters:

• element_coords (torch.Tensor \([N_e, N_b, D]\)) – Boundary element coordinates, where \(N_e\) = number of elements, \(N_b\) = number of basis functions, \(D\) = spatial dimension

• basis_order (int) – The order of the basis

• quadrature_order (int) – The order of quadrature rule to use

Returns:

For elements with uniform facets, a 5D tensor of shape \([N_e, N_f, N_q, D-1, D]\) where \(N_e\) is the number of elements, \(N_f\) the number of facets per element, \(N_q\) the number of quadrature points per facet, and \(D\) the spatial dimension.

For elements with mixed facets (Prism, Pyramid), a pair (tri_facet_jacobian, quad_facet_jacobian) of shapes \([N_e, N_{tf}, N_{tq}, D-1, D]\) and \([N_e, N_{qf}, N_{qq}, D-1, D]\), respectively.

Return type:

Tensor or Tuple[Tensor, Tensor]

classmethod element_to_facet(elements: Tensor, order: int) → Tensor | Tuple[Tensor, Tensor]

Maps element basis functions to facet basis functions.

For each facet of an element, extracts the subset of element basis functions that are active on that facet.

For mixed elements with both triangular and quadrilateral facets, returns separate tensors for each facet type.

For example:

# Hexahedral element (order=1)
hex_element = torch.randn(8)  # 8 basis functions
hex_facets = element_to_facet(hex_element, order=1)  # [6, 4] tensor

# Tetrahedral element (order=2)
tet_element = torch.randn(10)  # 10 basis functions
tet_facets = element_to_facet(tet_element, order=2)  # [4, 6] tensor

# Wedge element (order=2)
wedge_element = torch.randn(18)  # 18 basis functions
tri_facets, quad_facets = element_to_facet(wedge_element, order=2)
# tri_facets: [2, 6] tensor for triangular faces
# quad_facets: [3, 9] tensor for quadrilateral faces

Parameters:

• elements (Tensor) –

  2D tensor of shape \([N_e, N_b]\) containing element basis function values or 1D tensor of shape \([N_b]\) for a single element, where:

  • \(N_e\) = number of elements

  • \(N_b\) = number of basis functions

• order (int) – Order of the basis functions

Returns:

For elements with uniform facets (is_mix_facet is False), a tensor of shape \([N_e, N_f, N_b]\) (3D input) or \([N_f, N_b]\) (1D input), where \(N_e\) is the number of elements, \(N_f\) the number of facets, and \(N_b\) the number of basis functions per facet.

For elements with mixed facets, a pair (tri_facet, quad_facet) of shapes \([\ldots, N_{tf}, N_{tb}]\) and \([\ldots, N_{qf}, N_{qb}]\), respectively.

Return type:

Tensor or Tuple[Tensor, Tensor]

classmethod element_to_edge(elements: Tensor, order: int, unique: bool = True) → Tensor

Extract edge values from element values.

Maps element basis function values to edge basis function values using the edge connectivity.

Parameters:

• elements (Tensor) – ND tensor of shape […, N_b] containing basis function values, where N_b is the number of basis functions. The leading dimensions can be arbitrary (e.g., batch dimensions).

• order (int) – Order of the basis functions

• unique (bool, optional) – If True, returns only unique edges after sorting edge values. If False, returns all edges including duplicates. Default is True.

Returns:

edge – (N+1)D tensor of shape […, N_ed, N_eb] containing edge basis function values, where:

• […] = input tensor’s leading dimensions

• N_ed = number of edges

• N_eb = number of basis functions per edge

Return type:

Tensor

classmethod get_tri_mask() → Tensor

Returns a boolean mask indicating which faces are triangular.

For elements with mixed face types (triangular and quadrilateral faces), returns a boolean tensor marking which faces are triangular (True) vs quadrilateral (False).

The mask has length equal to the total number of faces and True values for each triangular face.

Examples

For a prism element with 2 triangular faces and 3 quadrilateral faces:

>>> Element.get_tri_mask()
tensor([True, True, False, False, False])

Returns:

tri_mask – 1D boolean Tensor of shape \([N_f]\) where \(N_f\) = number of facets

Return type:

Tensor

classmethod get_quad_mask() → Tensor

Returns a boolean mask indicating which faces are quadrilateral.

For elements with mixed face types (triangular and quadrilateral faces), returns a boolean tensor marking which faces are quadrilateral (True) vs triangular (False).

The mask has length equal to the total number of faces and True values for each quadrilateral face.

Examples

For a prism element with 2 triangular faces and 3 quadrilateral faces:

>>> Element.get_quad_mask()
tensor([False, False, True, True, True])

Returns:

quad_mask – 1D boolean Tensor of shape \([N_f]\) where \(N_f\) = number of facets

Return type:

Tensor

classmethod get_contour(order: int) → Tensor

Returns the contour points for a 2D element at the specified order.

For a 2D element, returns points along the element boundary at the given polynomial order. The points are ordered counterclockwise around the boundary.

Examples

Get contour points for a quadrilateral at order 2:

>>> quad = Quadrilateral()
>>> contour = quad.get_contour(2)
>>> contour.shape
torch.Size([8, 2])

Parameters:

order (int) – The polynomial order. Must be >= 1.

Returns:

contour – Tensor of shape \([N_p, 2]\) containing the \((x,y)\) coordinates of points along the element boundary, ordered counterclockwise, where \(N_p\) = number of points.

Return type:

Tensor

Notes

• Only implemented for 2D elements (dim=2)

• Points are distributed according to the polynomial order

• Points are ordered counterclockwise starting from the first vertex

The seven concrete shapes inherit the methods above; the subclass entries below mainly document their reference geometry (vertex coordinates, facet shape, type-string family).

class Line

Bases: Element

1D line reference element on \([0, 1]\).

The facet of every 2D element (edge of a triangle or quad) is a Line. Type strings: "line" (linear), "line3", "line4", …, "line11".

__init__()

class Triangle

Bases: Element

2D triangle reference element with vertices at \((0,0), (1,0), (0,1)\).

Simplex ordering: vertex nodes first, then edge nodes, then interior nodes. The facets are line segments. Type strings: "triangle" (linear), "triangle6", "triangle10", …, "triangle66".

__init__()

class Quadrilateral

Bases: Element

2D quadrilateral reference element on \([0, 1]^2\).

Tensor-product Lagrange nodes in lexicographic order — note this differs from the Gmsh/VTK boundary-spiral order, so loading external high-order quads must go through Element.reorder(). The facets are line segments. Type strings: "quad" (linear), "quad8", "quad9", "quad16", …, "quad121".

__init__()

class Tetrahedron

Bases: Element

3D tetrahedron reference element with vertices at \((0,0,0), (1,0,0), (0,1,0), (0,0,1)\).

Simplex ordering: vertex / edge / face / interior groupings, with indices within each group following the TensorMesh convention (see Gmsh / VTK ↔ TensorMesh node ordering). The facets are triangles. Type strings: "tetra" (linear), "tetra10", "tetra20", …, "tetra286".

__init__()

class Hexahedron

Bases: Element

3D hexahedron reference element on \([0, 1]^3\).

Tensor-product Lagrange nodes in lexicographic order — note this differs substantially from Gmsh/VTK, which walks along edges. The facets are quadrilaterals. Type strings: "hexahedron" (linear), "hexahedron20", "hexahedron27", "hexahedron64", …, "hexahedron1000".

__init__()

class Prism

Bases: Element

3D triangular prism (wedge) reference element — triangle base \(\times\) line height.

Mixed-facet element: two triangular caps + three quadrilateral sides (see Element.is_mix_facet). Note the string convention "wedge" differs from the class name Prism (carried over from meshio). Type strings: "wedge" (linear), "wedge18", "wedge40", "wedge75", …, "wedge550".

__init__()

class Pyramid

Bases: Element

3D pyramid reference element with a unit square base and apex at \((0, 0, 1)\).

Mixed-facet element: one quadrilateral base + four triangular sides (see Element.is_mix_facet). Type strings: "pyramid" (linear), "pyramid14".

__init__()

Type-string registry

Lookup tables and helpers that map between element type strings ("triangle", "tetra10", …) and the corresponding element class, spatial dimension, and polynomial order.

element_type2element(x: str) → Type[Element]

Resolve an element type string to its Element subclass.

The alphabetic prefix of x determines the shape ("triangle6" → Triangle, "hexahedron27" → Hexahedron, etc.); the trailing digits — which encode the basis-node count — are ignored here and instead consumed by element_type2order.

Parameters:

x (str) – Element type string. Must start with one of the supported prefixes: "line", "triangle", "quad", "tetra", "hexahedron", "pyramid", or "wedge".

Returns:

The Element subclass corresponding to x.

Return type:

Type[Element]

element_types: list[str]

List of every element type string the library understands — first-order shapes ("line", "triangle", "quad", "tetra", "hexahedron", "wedge", "pyramid") plus their higher-order counterparts ("triangle6", "quad9", "tetra10", "hexahedron27", "triangle10", …).

element_type2dimension: dict[str, int]

Map from element type string to spatial dimension ("line": 1, "triangle": 2, "tetra": 3, …).

element_type2order: dict[str, int]

Map from element type string to polynomial order ("triangle": 1, "triangle6": 2, "triangle10": 3, …).

Extension internals

These layers exist for people adding a new element type or a new quadrature rule. None of them are part of the stable user-facing API: their signatures may evolve between releases and their docstrings live in the source rather than here.

Polynomial space. tensormesh.element.Polynomial represents a single multivariate polynomial; tensormesh.element.Polynomials represents an arbitrarily-shaped batch. They are the building blocks of every shape-function set: get_basis_fns() solves a small Vandermonde system to obtain the Lagrange-interpolation coefficients on top of a chosen polynomial space (full poly_exp, tensor-product tens_exp, or the pyramid / prism ad-hoc spaces). To add a new element you typically override get_polynomial() to return the right polynomial space and let the base class produce the basis functions for you.

class Polynomial(exp: Tensor, coef: Tensor | None = None)

Bases: Module

A polynomial class representing multivariate polynomials of the form:

\[p(x_1,\ldots,x_n) = \sum_{i=1}^m c_i \prod_{j=1}^n x_j^{e_{ij}}\]

where:

• \(n\) is the number of variables (n_vars)

• \(m\) is the number of terms (n_terms)

• \(c_i\) are the coefficients (coef)

• \(e_{ij}\) are the exponents (exp)

• \(x_j\) are the variables

For example, the polynomial \(x^2y + 2xy^2 + 3\) has:

• n_vars = 2 (x, y)

• n_terms = 3

• coef = [1, 2, 3]

• exp = [[2, 1, 0], [1, 2, 0]] (rows index variables, columns index terms)

Examples

# Create polynomial with exponents and coefficients
exp = torch.tensor([[1, 2, 3], [2, 1, 0]]) # [n_vars(2), n_terms(3)]
coef = torch.tensor([1, 1, 1])             # [n_terms(3)]
poly = Polynomial(exp, coef)

# Print polynomial
print(poly)  # xy^2 + x^2y + x^3

# Evaluate polynomial at point x
x = torch.tensor([2, 3]) # [n_vars]
print(poly(x))  # 2*3^2 + 2^2*3 + 2^3 = 38

# Take derivatives
print(poly.deriv(0))  # y^2 + 2xy + 3x^2
print(poly.deriv(1))  # 2xy + x^2

__init__(exp: Tensor, coef: Tensor | None = None)

Initialize a polynomial with exponents and coefficients.

Parameters:

• exp (Tensor) –

  Exponents tensor of shape \([N_v, N_t]\) where:

  • \(N_v\) = number of variables

  • \(N_t\) = number of terms

• coef (Optional[Tensor], optional) – Coefficients tensor of shape \([N_t]\), by default None. If None, coefficients will be initialized as ones.

Examples

# Create polynomial x^2y + 2xy^2 + 3
exp = torch.tensor([[2,1,0], [1,2,0]])  # [2 vars, 3 terms]
coef = torch.tensor([1,2,3])            # [3 terms]
poly = Polynomial(exp, coef)

n_vars: int

Number of variables in polynomial :no-index:

n_terms: int

Number of terms in polynomial :no-index:

property device

property dtype

forward(x: Tensor) → Tensor

Evaluate the polynomial at given points.

Examples

# Create a polynomial p(x,y) = 2x + 3y^2
coef = torch.tensor([2.0, 3.0])
exp = torch.tensor([[1.0, 0.0], [0.0, 2.0]])
poly = Polynomial(exp, coef)

# Evaluate at single point
x = torch.tensor([1.0, 2.0])  # point (x=1, y=2)
poly(x)  # 2*1 + 3*2^2 = 14.0
# tensor(14.)

# Evaluate at multiple points
x = torch.tensor([[1.0, 2.0], [2.0, 1.0]])  # points (1,2) and (2,1)
poly(x)  # [2*1 + 3*2^2, 2*2 + 3*1^2] = [14.0, 7.0]
# tensor([14., 7.])

Parameters:

x (Tensor) – Input tensor of shape [n_batch, n_vars] or [n_vars].

Returns:

Evaluated polynomial values of shape [n_batch] or scalar.

Return type:

Tensor

get_exp_terms(x: Tensor) → Tensor

Compute exponential terms for polynomial evaluation.

For a polynomial with terms \(c_i x_1^{e_{i1}} x_2^{e_{i2}} \cdots x_n^{e_{in}}\), computes the exponential terms \(x_1^{e_{i1}} x_2^{e_{i2}} \cdots x_n^{e_{in}}\) for each term i.

Examples

# Create polynomial p(x,y) = 2x^2y + 3xy^2
exp = torch.tensor([[2,1], [1,2]])  # exponents for each term
coef = torch.tensor([2.0, 3.0])     # coefficients
poly = Polynomial(exp, coef)

# Single point evaluation
x = torch.tensor([2.0, 3.0])  # point (x=2, y=3)
terms = poly.get_exp_terms(x)  # [2^2 * 3^1, 2^1 * 3^2]
# tensor([12., 54.])

# Multiple point evaluation
x = torch.tensor([[2.0, 3.0], [1.0, 2.0]])  # points (2,3) and (1,2)
terms = poly.get_exp_terms(x)  # [[2^2 * 3^1, 2^1 * 3^2], [1^2 * 2^1, 1^1 * 2^2]]
# tensor([[12., 54.],
#         [ 2.,  4.]])

Parameters:

x (Tensor) – Input points tensor of shape: * [n_vars] for single point * [n_batch, n_vars] for multiple points where: * n_vars = number of variables * n_batch = number of points

Returns:

Exponential terms tensor of shape: * [n_terms] for single point input * [n_batch, n_terms] for multiple point input where: * n_terms = number of polynomial terms

Return type:

Tensor

deriv(var_ind: int = 0) → Polynomial

Compute the derivative of the polynomial with respect to a variable.

For a polynomial \(p(x,y,z) = ax^ny^mz^k\), the derivative with respect to x is: \(\frac{\partial p}{\partial x} = nax^{n-1}y^mz^k\) if n>0, or 0 if n=0

Examples

# Create polynomial p(x,y) = x^2y + 2xy^2 + 3
exp = torch.tensor([[2, 1, 0], [1, 2, 0]]) # [n_vars(2), n_terms(3)]
coef = torch.tensor([1, 2, 3])             # [n_terms(3)]
poly = Polynomial(exp, coef)

# Take derivative with respect to x
dx = poly.deriv(0)  # 2xy + 2y^2
print(dx)

# Take derivative with respect to y
dy = poly.deriv(1)  # x^2 + 4xy
print(dy)

# Evaluate derivatives at point (2,3)
x = torch.tensor([2.0, 3.0])
print(dx(x))  # 2*2*3 + 2*3^2 = 30
print(dy(x))  # 2^2 + 4*2*3 = 28

Parameters:

var_ind (int) – the index of the variable to be differentiated

Returns:

deriv_poly – the derivative of the polynomial

Return type:

Polynomial

grad() → Polynomials

Compute the gradient of the polynomial.

For a polynomial \(p(x_1,\ldots,x_n)\), returns a vector of partial derivatives: \(\nabla p = [\frac{\partial p}{\partial x_1}, \ldots, \frac{\partial p}{\partial x_n}]\)

For example, given \(p(x,y) = ax^ny^m\), the gradient is: \(\nabla p = [nax^{n-1}y^m, max^ny^{m-1}]\)

Examples

# Create polynomial p(x,y) = x^2y + 2xy^2 + 3
exp = torch.tensor([[2, 1, 0], [1, 2, 0]]) # [n_vars(2), n_terms(3)]
coef = torch.tensor([1, 2, 3])             # [n_terms(3)]
poly = Polynomial(exp, coef)

# Take gradient
grad = poly.grad()  # [2xy + 2y^2, x^2 + 4xy]
print(grad)

# Evaluate gradient at point (2,3)
x = torch.tensor([2.0, 3.0])
print(grad(x))  # [30, 28]

Parameters:

None

Returns:

grad_poly – the gradient of the polynomial

Return type:

Polynomials

reset_coef(coef: Tensor)

Reset the coefficients of the polynomial while keeping the exponents unchanged.

Parameters:

coef (Tensor) – The new coefficients to set. Must match the shape, device and dtype of the existing coefficients.

Returns:

self – Returns self for method chaining.

Return type:

Polynomial

repeat(*args) → Polynomials

Repeat the polynomial along specified dimensions.

Creates a new Polynomials object by repeating the current polynomial’s coefficients and exponents according to the provided repeat dimensions. This is similar to torch.repeat().

For example:

• repeat(3) creates 3 copies along a new first dimension

• repeat(2,3) creates a 2x3 grid of copies in new dimensions

The coefficients and exponents are repeated while preserving the polynomial structure, effectively creating multiple independent copies of the same polynomial.

Examples

# Create a polynomial
exp = torch.tensor([[2, 1], [1, 0]])  # x^2y + x
coef = torch.tensor([1.0, 1.0])
poly = Polynomial(exp, coef)

# Repeat 3 times along new first dimension
polys1 = poly.repeat(3)  # Shape: [3, n_vars, n_terms]

# Create 2x3 grid of copies
polys2 = poly.repeat(2, 3)  # Shape: [2, 3, n_vars, n_terms]

# Evaluate repeated polynomials
x = torch.randn(2, 3, 2)  # [2, 3, n_vars]
y = polys2(x)  # [2, 3] outputs

Parameters:

*args (int) – The number of repetitions for each new dimension. Similar to torch.repeat() arguments.

Returns:

A new Polynomials object with the repeated structure. The shape will be [*args, *original_shape].

Return type:

Polynomials

classmethod lin_exp(n_vars: int, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomial

Creates a linear polynomial with n_vars variables. Creates a polynomial with n_vars+1 terms:

• A constant term (all exponents 0)

• Linear terms for each variable (exponent 1 for that variable, 0 for others)

For example, with n_vars=2, creates polynomial with terms: \(1, x, y\)

Examples

# Create linear polynomial with 2 variables (1 + x + y)
poly = Polynomial.lin_exp(2)
print(poly)  # 1 + x + y

# Create linear polynomial with 3 variables (1 + x + y + z)
poly = Polynomial.lin_exp(3)
print(poly)  # 1 + x + y + z

# Specify dtype and device
poly = Polynomial.lin_exp(2, dtype=torch.float64, device=torch.device('cuda'))
print(poly.dtype)  # torch.float64
print(poly.device)  # cuda:0

Parameters:

• n_vars (int) – number of vairables

• dtype (dtype) – the data type of the polynomial

• device (device) – the device of the polynomial

Returns:

polynomial

Return type:

Polynomial n_vars=n_vars n_terms=n_vars+1

classmethod poly_exp(n_vars: int, dim: int, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomial

Creates a polynomial with terms up to a maximum total degree.

For n_vars variables and maximum degree dim, includes all terms where the sum of exponents is less than or equal to dim, sorted by total degree.

For example, with n_vars=2, dim=2, creates polynomial with terms: \(1, x, y, x^2, xy, y^2\)

For n_vars=1, creates polynomial with terms up to power dim: \(1, x, x^2, \ldots, x^{dim}\)

Examples

# Create polynomial with terms up to degree 2 in 2 variables
poly = Polynomial.poly_exp(2, 2)
print(poly)  # 1 + x + y + x^2 + xy + y^2

# Create polynomial with terms up to degree 3 in 1 variable
poly = Polynomial.poly_exp(1, 3)
print(poly)  # 1 + x + x^2 + x^3

# Create polynomial with terms up to degree 2 in 3 variables
poly = Polynomial.poly_exp(3, 2)
print(poly)  # 1 + x + y + z + x^2 + xy + xz + y^2 + yz + z^2

# Specify dtype and device
poly = Polynomial.poly_exp(2, 2, dtype=torch.float64, device=torch.device('cuda'))
print(poly.dtype)  # torch.float64
print(poly.device)  # cuda:0

Parameters:

• n_vars (int) – number of vairables

• dim (int) – number of maximum dimension

• dtype (dtype) – the data type of the polynomial

• device (device) – the device of the polynomial

Returns:

polynomial

Return type:

Polynomial n_vars=n_vars

classmethod tens_exp(n_vars: int, dim: int, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomial

Creates a tensor product polynomial with terms up to maximum degree in each variable.

For n_vars variables and maximum degree dim, includes all terms where each exponent is less than or equal to dim, sorted by total degree.

For example, with n_vars=2, dim=2, creates polynomial with terms: \(P(x,y) = a_0 + a_1x + a_2y + a_3x^2 + a_4xy + a_5y^2 + a_6x^2y + a_7xy^2 + a_8x^2y^2\)

For n_vars=1, dim=d, creates polynomial with terms up to power dim: \(P(x) = a_0 + a_1x + a_2x^2 + ... + a_dx^d\)

Examples

# Create tensor product polynomial with 2 variables up to degree 2
poly = Polynomial.tens_exp(2, 2)
print(poly)  # 1 + x + y + x^2 + xy + y^2 + x^2y + xy^2 + x^2y^2

# Create tensor product polynomial with 1 variable up to degree 3
poly = Polynomial.tens_exp(1, 3)
print(poly)  # 1 + x + x^2 + x^3

# Specify dtype and device
poly = Polynomial.tens_exp(2, 2, dtype=torch.float64, device=torch.device('cuda'))
print(poly.dtype)  # torch.float64
print(poly.device)  # cuda:0

Parameters:

• n_vars (int) – number of vairables

• dim (int) – number of maximum dimension

• dtype (dtype) – the data type of the polynomial

• device (device) – the device of the polynomial

Returns:

polynomial

Return type:

Polynomial n_vars=n_vars n_terms=(dim+1)**n_vars

classmethod pyr_exp(order: int, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomial

Creates a pyramid polynomial with terms up to maximum total degree.

For order n, includes all terms where x+z < n+1 and y+z < n+1, sorted by total degree.

For example, with order=1, creates polynomial with terms: \(P(x,y,z) = a_0 + a_1x + a_2y + a_3z\)

With order=2, creates polynomial with terms: \(P(x,y,z) = a_0 + a_1x + a_2y + a_3z + a_4x^2 + a_5xy + a_6y^2 + a_7xz + a_8yz + a_9z^2\)

Examples

# Create order 1 pyramid polynomial
poly = Polynomial.pyr_exp(1)
print(poly)  # 1 + x + y + z

# Create order 2 pyramid polynomial with custom dtype and device
poly = Polynomial.pyr_exp(2, dtype=torch.float64, device=torch.device('cuda'))
print(poly.dtype)  # torch.float64
print(poly.device)  # cuda:0

# Evaluate polynomial at point x
x = torch.tensor([1.0, 2.0, 0.5], device='cuda', dtype=torch.float64)  # [x,y,z]
print(poly(x))  # Evaluates polynomial at (1,2,0.5)

Parameters:

• order (int) – Maximum polynomial order

• dtype (dtype, optional) – Data type of polynomial coefficients

• device (device, optional) – Device to store polynomial on

Returns:

polynomial – Pyramid polynomial with n_vars=3 and appropriate number of terms

Return type:

Polynomial

classmethod pri_exp(order: int, dtype: dtype = torch.float32, device: device = torch.device('cpu')) → Polynomial

Creates a prismatic polynomial with terms up to maximum total degree.

For order n, includes all terms where x+y < n+1, sorted by total degree.

For example, with order=1, creates polynomial with terms: \(P(x,y,z) = a_0 + a_1x + a_2y + a_3z\)

With order=2, creates polynomial with terms: \(P(x,y,z) = a_0 + a_1x + a_2y + a_3z + a_4x^2 + a_5xy + a_6y^2 + a_7xz + a_8yz + a_9z^2\)

Examples

# Create prismatic polynomial of order 2
poly = Polynomial.pri_exp(2)
print(poly)  # 1 + x + y + z + x^2 + xy + y^2 + xz + yz + z^2

# Create with specific dtype and device
poly = Polynomial.pri_exp(2, dtype=torch.float64, device=torch.device('cuda'))
print(poly.dtype)  # torch.float64
print(poly.device)  # cuda:0

# Evaluate polynomial at point x
x = torch.tensor([1.0, 2.0, 0.5], device='cuda', dtype=torch.float64)  # [x,y,z]
print(poly(x))  # Evaluates polynomial at (1,2,0.5)

Parameters:

• order (int) – Maximum polynomial order

• dtype (dtype, optional) – Data type of polynomial coefficients

• device (device, optional) – Device to store polynomial on

Returns:

polynomial – Prismatic polynomial with n_vars=3 and appropriate number of terms

Return type:

Polynomial

class Polynomials(exp: Tensor, coef: Tensor | None = None)

Bases: Module

A collection of polynomials that can be evaluated and manipulated together.

This class represents a batch of polynomials, allowing vectorized operations across multiple polynomials. Each polynomial has the same number of variables and terms, but can have different coefficients.

The polynomials are stored in a batched format with coefficients and exponents tensors:

• Coefficients tensor shape: [n_poly1, …, n_terms]

• Exponents tensor shape: [n_poly1, …, n_vars, n_terms]

Where n_poly1, … are the batch dimensions.

Examples

# Create a batch of 2 quadratic polynomials in 2 variables
exp = torch.tensor([[[0,1,0,2,0], [0,0,1,0,2]],
                  [[0,1,0,2,0], [0,0,1,0,2]]])  # [2, 2, 5]
coef = torch.tensor([[1,2,3,4,5], [2,3,4,5,6]]) # [2, 5]
polys = Polynomials(exp, coef)

# Print batch shape and dimensions
print(polys.shape)     # (2,)
print(polys.n_vars)    # 2
print(polys.n_terms)   # 5

# Evaluate polynomials at points
x = torch.tensor([[1.0, 2.0], [3.0, 4.0]])  # [2, 2]
y = polys(x)  # Evaluates both polynomials at their respective points

# Take derivatives
d_polys = polys.deriv(0)  # Derivative with respect to first variable

# Convert to individual polynomials
poly_list = [polys[i] for i in range(len(polys))]

# Create with specific dtype and device
exp = exp.cuda().double()
coef = coef.cuda().double()
polys = Polynomials(exp, coef)
print(polys.dtype)    # torch.float64
print(polys.device)   # cuda:0

Parameters:

• exp (Tensor) – Tensor of exponents with shape [n_poly1, …, n_vars, n_terms]

• coef (Tensor, optional) – Tensor of coefficients with shape [n_poly1, …, n_terms]. If None, defaults to ones.

_coef

The coefficients tensor

Type:

Tensor

_exp

The exponents tensor

Type:

Tensor

n_polys

The batch dimensions

Type:

Tuple[int, …]

n_vars

Number of variables in each polynomial

Type:

int

n_terms

Number of terms in each polynomial

Type:

int

device

The device the tensors are stored on

Type:

device

dtype

The data type of the tensors

Type:

dtype

shape

The batch dimensions (same as n_polys)

Type:

Tuple[int, …]

__init__(exp: Tensor, coef: Tensor | None = None)

Initialize multiple polynomials with exponents and coefficients.

Parameters:

• exp (Tensor) – Exponents tensor of shape \([N_1, \ldots, N_v, N_t]\) where: * \(N_1, \ldots\) = batch dimensions for multiple polynomials * \(N_v\) = number of variables in each polynomial * \(N_t\) = number of terms in each polynomial

• coef (Optional[Tensor], optional) – Coefficients tensor of shape \([N_1, \ldots, N_t]\) where: * \(N_1, \ldots\) = batch dimensions for multiple polynomials * \(N_t\) = number of terms in each polynomial If None, coefficients will be initialized as ones.

Examples

# Create 2x3 grid of polynomials x^2y + 2xy^2 + 3
exp = torch.tensor([[[2,1,0], [1,2,0]]]).repeat(2,3,1,1)  # [2,3,2,3]
coef = torch.tensor([[1,2,3]]).repeat(2,3,1)              # [2,3,3]
polys = Polynomials(exp, coef)  # 2x3 grid of polynomials

# Print shape
print(polys.shape)  # (2,3)

# Access individual polynomials
print(polys[0,0])  # x^2y + 2xy^2 + 3
print(polys[1,2])  # x^2y + 2xy^2 + 3

numel()

property ndim

dim()

forward(x: Tensor) → Tensor

Evaluates the polynomial at the given input points.

For a polynomial \(p(x) = \sum_i c_i \prod_j x_j^{e_{ij}}\), computes:

1. Evaluates each term by raising inputs to exponents

2. Multiplies terms by coefficients

3. Sums the terms

For example, for polynomial \(2x^2y + 3xy^2\) with point \([2,3]\):

\[2(2^2 \cdot 3^1) + 3(2^1 \cdot 3^2) = 2(12) + 3(18) = 78\]

The computation is vectorized across batches and/or multiple polynomials.

Examples

# Single polynomial evaluation
exp = torch.tensor([[2, 1], [1, 2]])  # x^2y + xy^2
coef = torch.tensor([2, 3])           # 2x^2y + 3xy^2
poly = Polynomial(exp, coef)
x = torch.tensor([2.0, 3.0])          # Point [2,3]
print(poly(x))                        # 2*(2^2*3) + 3*(2*3^2) = 78.0

# Batch evaluation
x_batch = torch.tensor([[2.0, 3.0], [1.0, 2.0]])  # 2 points
print(poly(x_batch))                  # [78.0, 14.0]

# Multiple polynomials
polys = Polynomials.stack([poly, poly])     # 2 copies of same polynomial
print(polys(x))                       # [78.0, 78.0]

# Batch + multiple polynomials
print(polys(x_batch))                 # [[78.0, 78.0], [14.0, 14.0]]

Parameters:

x (Tensor) – Input points of shape [n_batch, n_poly1, ..., n_vars] or [n_poly1, ..., n_vars].

Returns:

Evaluated values of shape [n_batch, n_poly1, ...] or [n_poly1, ...].

Return type:

Tensor

get_exp_terms(x: Tensor) → Tensor

Evaluates the polynomial terms for each input point by raising to the appropriate exponents.

For each input point \(x\) and term with exponents \(e\), computes:

\[x_0^{e_0} \cdot x_1^{e_1} \cdot \ldots \cdot x_n^{e_n}\]

For example, for polynomial \(x^2y\) with point \([2,3]\), computes:

\[2^2 \cdot 3^1 = 12\]

The computation is vectorized across batches and/or multiple polynomials.

Examples

# Single polynomial evaluation
exp = torch.tensor([[2, 1], [1, 2]])  # x^2y + xy^2
coef = torch.tensor([2, 3])           # 2x^2y + 3xy^2
poly = Polynomial(exp, coef)
x = torch.tensor([2.0, 3.0])          # Point [2,3]
terms = poly.get_exp_terms(x)         # [12, 18] (2^2*3, 2*3^2)

# Batch evaluation
x_batch = torch.tensor([[2.0, 3.0], [1.0, 2.0]])  # 2 points
terms = poly.get_exp_terms(x_batch)   # [[12, 18], [2, 4]]

# Multiple polynomials
polys = Polynomials.stack([poly, poly])     # 2 copies of same polynomial
terms = polys.get_exp_terms(x)        # [[12, 18], [12, 18]]

# Batch + multiple polynomials
terms = polys.get_exp_terms(x_batch)  # [[[12, 18], [12, 18]],
                                     #  [[2, 4], [2, 4]]]

Parameters:

x (Tensor) – Input points of shape [n_batch, n_poly1, ..., n_vars] or [n_poly1, ..., n_vars]. Can include a batch dimension. n_poly1, ... are optional polynomial batch dimensions. Last dimension must match number of variables.

Returns:

Output of shape [n_batch, n_poly1, ..., n_terms] or [n_poly1, ..., n_terms]. Evaluated terms for each input point. Output has same batch/polynomial dimensions as input. Last dimension contains values for each term.

Return type:

Tensor

apply_coefficient(x: Tensor) → Tensor

Applies the polynomial coefficients to the evaluated terms.

For each polynomial with coefficients \(c\) and evaluated terms \(t\), computes:

\[\sum_i c_i t_i\]

For example, for polynomial \(2x^2 + 3y\) with evaluated terms \([4,3]\), computes:

\[2 \cdot 4 + 3 \cdot 3 = 17\]

The computation is vectorized across batches and/or multiple polynomials.

The coefficients are broadcast to match the batch dimensions of the input.

Examples

# Single polynomial evaluation
exp = torch.tensor([[2, 1], [0, 1]])  # x^2, y
coef = torch.tensor([2, 3])           # 2x^2 + 3y
poly = Polynomial(exp, coef)
terms = torch.tensor([4, 3])          # terms = [x^2=4, y=3]
poly.apply_coefficient(terms)          # 2*4 + 3*3 = 17

# Batch evaluation
terms = torch.tensor([[4, 3], [1, 2]]) # 2 points
poly.apply_coefficient(terms)          # [2*4 + 3*3, 2*1 + 3*2]

# Multiple polynomials
exp = torch.tensor([[[2, 1], [0, 1]], # 2x^2 + 3y
                  [[1, 2], [1, 0]]])  # x^2y + xy
coef = torch.tensor([[2, 3], [1, 1]])
poly = Polynomials(exp, coef)
terms = torch.tensor([[4, 3], [12, 2]]) # terms for each polynomial
poly.apply_coefficient(terms)          # [2*4 + 3*3, 1*12 + 1*2]

Parameters:

x (Tensor) – ND Tensor of shape [n_batch, n_poly1, …, n_terms] or [n_poly1, …, n_terms]

Returns:

ND Tensor of shape [n_batch, n_poly1, …, n_poly2] or [n_poly1, …, n_poly2] where:

• n_batch = number of input points (optional)

• n_poly1, …, n_poly2 = polynomial dimensions

Return type:

Tensor

map(x: Tensor) → Tensor

Evaluates the polynomial at given input points.

For each polynomial with coefficients c and exponents e, computes:

\[\sum_i c_i \prod_j x_j^{e_{ij}}\]

For example, for polynomial \(2x^2 + 3y\) evaluated at point \((2,3)\), computes:

\[2 \cdot 2^2 + 3 \cdot 3 = 17\]

The computation is vectorized across batches and/or multiple polynomials. The input points are broadcast to match the polynomial dimensions.

Examples

# Single polynomial evaluation
exp = torch.tensor([[2, 1], [1, 2]])  # x^2y, xy^2
coef = torch.tensor([2, 3])           # 2x^2y + 3xy^2
poly = Polynomial(exp, coef)
x = torch.tensor([2.0, 3.0])          # Point (2,3)
poly.map(x)                           # 2*(2^2*3) + 3*(2*3^2) = 24 + 54 = 78
# tensor(78.)

# Batch evaluation
x = torch.tensor([[2.0, 3.0], [1.0, 2.0]])  # Two points
poly.map(x)                                  # Evaluate at both points
# tensor([78., 14.])

# Multiple polynomials
polys = Polynomials.stack([poly, poly])  # Two copies of same polynomial
x = torch.tensor([2.0, 3.0])       # Single point
polys.map(x)                       # Evaluate both polynomials
# tensor([78., 78.])

# Batch + multiple polynomials
x = torch.tensor([[2.0, 3.0], [1.0, 2.0]])  # Two points
polys.map(x)                                 # Evaluate both polys at both points
# tensor([[78., 78.],
#         [14., 14.]])

Parameters:

x (Tensor) – Input points tensor of shape [n_batch, n_vars] or [n_vars]. Each row represents a point to evaluate the polynomial(s) at.

Returns:

Evaluated polynomial values. Shape is: - [n_batch, n_poly1, …] if input is [n_batch, n_vars] - [n_poly1, …] if input is [n_vars] Where n_poly1, … are the polynomial batch dimensions.

Return type:

Tensor

map_exp_terms(x: Tensor) → Tensor

Evaluates the polynomial terms at given input points by applying exponents.

For each polynomial with exponents e, computes:

\[\prod_j x_j^{e_{ij}}\]

For example, for polynomial terms \(x^2y, xy^2\) evaluated at point \((2,3)\), computes:

\[[2^2 \cdot 3^1, 2^1 \cdot 3^2] = [12, 18]\]

The computation is vectorized across batches and/or multiple polynomials. The input points are broadcast to match the polynomial dimensions.

Examples

# Single polynomial, single point
exp = torch.tensor([[2, 1], [1, 2]])  # x^2y, xy^2
poly = Polynomial(exp)
x = torch.tensor([2.0, 3.0])  # Point (2,3)
terms = poly.map_exp_terms(x)  # Evaluate terms
# tensor([12., 18.])  # 2^2 * 3^1, 2^1 * 3^2

# Single polynomial, batch of points
x = torch.tensor([[2.0, 3.0], [1.0, 2.0]])  # Two points
terms = poly.map_exp_terms(x)  # Evaluate terms at each point
# tensor([[12., 18.],  # Point 1: 2^2 * 3^1, 2^1 * 3^2
#         [1., 4.]])   # Point 2: 1^2 * 2^1, 1^1 * 2^2

# Multiple polynomials
polys = Polynomials.stack([poly, poly])  # Two copies
x = torch.tensor([2.0, 3.0])       # Single point
terms = polys.map_exp_terms(x)      # Evaluate terms for both polys
# tensor([[12., 18.],
#         [12., 18.]])

# Batch + multiple polynomials
x = torch.tensor([[2.0, 3.0], [1.0, 2.0]])  # Two points
terms = polys.map_exp_terms(x)               # Evaluate at each point
# tensor([[[12., 18.],
#          [12., 18.]],
#         [[1., 4.],
#          [1., 4.]]])

Parameters:

x (Tensor) – Input points tensor of shape [n_batch, n_vars] or [n_vars]. Each row represents a point to evaluate the polynomial terms at.

Returns:

Evaluated polynomial term values. Shape is: - [n_batch, n_poly1, …, n_terms] if input is [n_batch, n_vars] - [n_poly1, …, n_terms] if input is [n_vars] Where n_poly1, … are the polynomial batch dimensions.

Return type:

Tensor

reshape(*args: Sequence[int]) → Polynomials

Reshapes the polynomial batch dimensions.

Creates a new Polynomials object with the coefficients and exponents reshaped to the specified dimensions. The total number of elements must remain the same.

Similar to torch.reshape(), this operation changes the batch dimensions while preserving the polynomial structure. The n_vars and n_terms dimensions are preserved at the end.

Parameters:

*args (Sequence[int]) – The new shape dimensions. The product of these dimensions must equal the product of the original batch dimensions.

Returns:

A new Polynomials object with reshaped batch dimensions.

Return type:

Polynomials

Examples

>>> poly = Polynomials(exp, coef)  # shape (6,)
>>> reshaped = poly.reshape(2,3)   # shape (2,3)

transpose(*args: Sequence[int]) → Polynomials

Transposes the polynomial batch dimensions.

Creates a new Polynomials object with the coefficients and exponents transposed according to the specified dimension ordering. The n_vars and n_terms dimensions are preserved at the end.

Similar to torch.transpose(), this operation permutes the batch dimensions while preserving the polynomial structure.

Parameters:

*args (Sequence[int]) – The new ordering of dimensions. Must include all dimensions up to dim().

Returns:

A new Polynomials object with transposed batch dimensions.

Return type:

Polynomials

Examples

>>> poly = Polynomials(exp, coef)  # shape (2,3)
>>> transposed = poly.transpose(1,0)  # shape (3,2)

deriv(var_ind: int = 0) → Polynomials

Compute the derivative of the polynomial with respect to a variable.

For a polynomial \(p(x_1,\ldots,x_n)\), computes \(\frac{\partial p}{\partial x_i}\) where i is the specified variable index.

For example, given \(p(x,y) = ax^ny^m\), the derivatives are:

• \(\frac{\partial p}{\partial x} = nax^{n-1}y^m\)

• \(\frac{\partial p}{\partial y} = max^ny^{m-1}\)

The derivative is computed by:

1. Decrementing the exponent of the specified variable by 1

2. Multiplying coefficients by the original exponent

3. Setting terms with exponent 0 to 0

Examples

# Single polynomial derivative
exp = torch.tensor([[2, 1], [1, 2]])  # x^2y + xy^2
coef = torch.tensor([2, 3])           # 2x^2y + 3xy^2
poly = Polynomial(exp, coef)

# Derivative with respect to x
dx = poly.deriv(0)                    # 4xy + 3y^2

# Derivative with respect to y
dy = poly.deriv(1)                    # 2x^2 + 6xy

# Multiple polynomials
polys = Polynomials.stack([poly, poly])     # [2x^2y + 3xy^2, 2x^2y + 3xy^2]
dx = polys.deriv(0)                   # [4xy + 3y^2, 4xy + 3y^2]
dy = polys.deriv(1)                   # [2x^2 + 6xy, 2x^2 + 6xy]

Parameters:

var_ind (int) – the index of the variable to be differentiated

Returns:

deriv_poly – the derivative of the polynomial

Return type:

Polynomial

grad() → Polynomials

Compute the gradient of the polynomial.

For a polynomial \(p(x_1,\ldots,x_n)\), returns a vector of partial derivatives: \(\nabla p = [\frac{\partial p}{\partial x_1}, \ldots, \frac{\partial p}{\partial x_n}]\)

For example, given \(p(x,y) = ax^ny^m\), the gradient is: \(\nabla p = [nax^{n-1}y^m, max^ny^{m-1}]\)

The gradient is computed by taking the derivative with respect to each variable:

1. For each variable i=1…n: - Compute \(\frac{\partial p}{\partial x_i}\) using deriv(i)

2. Stack the derivatives into a vector

Examples

# Single polynomial gradient
exp = torch.tensor([[2, 1], [1, 2]])  # x^2y + xy^2
coef = torch.tensor([2, 3])           # 2x^2y + 3xy^2
poly = Polynomial(exp, coef)
grad = poly.grad()                    # [4xy + 3y^2, 2x^2 + 6xy]

# Multiple polynomials
polys = Polynomials.stack([poly, poly])     # [2x^2y + 3xy^2, 2x^2y + 3xy^2]
grad = polys.grad()                   # [[4xy + 3y^2, 4xy + 3y^2],
                                     #  [2x^2 + 6xy, 2x^2 + 6xy]]

# Evaluate gradient at points
x = torch.tensor([2.0, 3.0])          # Point [2,3]
grad_vals = grad(x)                   # [[36, 36], [24, 24]]

# Batch evaluation
x_batch = torch.tensor([[2.0, 3.0],   # 2 points
                      [1.0, 2.0]])
grad_vals = grad(x_batch)             # [[[36, 36], [24, 24]],
                                     #  [[8, 8], [4, 4]]]

Parameters:

None

Returns:

The gradient of the polynomial, with shape [n_vars, n_poly1, ..., n_polyn], n_vars variables and n_terms terms each.

Return type:

Polynomials

reset_coef(coef: Tensor)

Reset the coefficients of the polynomial

Parameters:

coef (Tensor) – Coefficients of shape [n_poly, n_terms].

repeat(*args)

Repeat the polynomial along specified dimensions.

Creates a new Polynomials object by repeating the current polynomial’s coefficients and exponents according to the provided repeat dimensions. This is similar to torch.repeat().

For example: - repeat(3) creates 3 copies along a new first dimension - repeat(2,3) creates a 2x3 grid of copies in new dimensions

The coefficients and exponents are repeated while preserving the polynomial structure, effectively creating multiple independent copies of the same polynomial.

Parameters:

*args (int) – The number of repetitions for each dimension. Must match the number of dimensions in the polynomial (self.dim()).

Returns:

A new Polynomials object with the repeated structure. The shape will be [args[0]*original_shape[0], args[1]*original_shape[1], …].

Return type:

Polynomials

classmethod stack(polys: Sequence[Polynomial | Polynomials]) → Polynomials

Stack multiple polynomials into a single Polynomials object.

Takes a sequence of Polynomial or Polynomials objects and stacks them along a new first dimension. All polynomials must have the same number of variables and terms.

Examples

# Create two simple polynomials: x^2 + y and 2x + y^2
exp1 = torch.tensor([[2,1], [0,1]])  # Exponents for x^2 + y
coef1 = torch.tensor([1.0, 1.0])     # Coefficients [1,1]
p1 = Polynomial(exp1, coef1)

exp2 = torch.tensor([[1,0], [0,2]])  # Exponents for 2x + y^2
coef2 = torch.tensor([2.0, 1.0])     # Coefficients [2,1]
p2 = Polynomial(exp2, coef2)

# Stack the polynomials
stacked = Polynomials.stack([p1, p2])
print(stacked.shape)  # (2,)

# Evaluate stacked polynomials at point [2,3]
x = torch.tensor([2.0, 3.0])
print(stacked(x))  # [7, 13] = [(2^2 + 3), (2*2 + 3^2)]

Parameters:

polys (Sequence[Union[Polynomial, Polynomials]]) – Sequence of polynomials to stack. All must have same n_vars, n_terms, dtype and device.

Returns:

A new Polynomials object with shape [len(polys), …] containing the stacked polynomials.

Return type:

Polynomials

Examples

>>> p1 = Polynomial(exp1, coef1)  # First polynomial
>>> p2 = Polynomial(exp2, coef2)  # Second polynomial with same structure
>>> stacked = Polynomials.stack([p1, p2])
>>> stacked.shape
(2,)

Basis-node and quadrature factories. The reference-element interpolation-node coordinates and quadrature rules live in three internal files; read the source if you need to understand exactly how get_basis() or get_quadrature() are implemented for each shape:

• tensormesh/element/basis.py — lin_basis, tri_basis, quad_basis, tet_basis, hex_basis, pyr_basis, pri_basis and the facet-basis index helpers.

• tensormesh/element/quadrature.py — lin_quadrature, tri_quadrature, quad_quadrature, tet_quadrature, hex_quadrature, pyr_quadrature, pri_quadrature plus the facet-quadrature variants.

• tensormesh/element/normal.py — outwards_normal_2d / outwards_normal_3d, used internally to compute the outward facet normals consumed by nanson_scale.

The supported extension path is to subclass Element and override its basis / quadrature / facet hooks. Use the existing seven shapes as templates — each subclass is short and self-contained.