Hvfem module (hvfem.h)

Functions

PetscErrorCode computeCellJacobian(Cell *cell)

Computes a cell’s Jacobian matrix, its inverse, and determinant.

The function computes the Jacobian of the affine transformation that maps points from the reference tetrahedron (vertices at (0,0,0), (1,0,0), (0,1,0), (0,0,1)) to the physical tetrahedron defined by cell->coordinates.

Steps performed:

Constructs the Jacobian matrix using differences of vertex coordinates.
Computes the determinant of the Jacobian.
Computes the cofactor matrix.
Computes the adjugate matrix (transpose of the cofactor matrix).
Computes the inverse Jacobian by dividing the adjugate by the determinant.

This function assumes a non-degenerate tetrahedron (non-zero volume). The Jacobian can be used for transforming gradients and for integration over the physical tetrahedron.

Parameters:

cell – [inout] Cell whose jacobian/invJacobian/detJacobian are filled.
cell – [inout] Pointer to a Cell structure containing:
- coordinates: array of size 12 with the tetrahedron vertex coordinates in the order [x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3].
- jacobian: 3x3 matrix to store the computed Jacobian.
- invJacobian: 3x3 matrix to store the inverse Jacobian.
- detJacobian: scalar to store the determinant.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS always.

PetscErrorCode computeCellOrientation(Cell *cell)

Computes the orientation of faces and edges for a tetrahedral cell.

The function extracts the orientation of each face and edge based on the cell’s transitive closure. The convention for indexing is as follows:

Faces: indices start at position 2 in closure, orientations at position 3, total of NUM_FACES_PER_CELL faces.
Edges: indices start after faces (2 + NUM_FACES_PER_CELL*2), orientations at the next position, total of NUM_EDGES_PER_CELL edges.

The extracted orientation values are then mapped to the PETGEM convention for basis functions:

Faces: [-3,-2,-1,0,1,2] → [4,3,5,0,1,2]
Edges: negative → -1, positive → 1

This allows the shape function construction routines to correctly handle local-to-global orientation transformations.

Parameters:

cell – [inout] Cell whose orientation (faces, edgeSigns) is filled.
cell – [inout] Pointer to a Cell structure containing:
- closure: array representing the cell’s transitive closure (e.g., from vertices to edges/faces).
- orientation: array to store the computed face and edge orientations.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS always.

PetscErrorCode computeNum1DQuadraturePoints(const PetscInt nord, Quadrature1D *quadrature)

Determines the number of 1D Gauss-Legendre quadrature points.

Computes the minimum number of 1D Gauss-Legendre points required to integrate polynomials up to degree 2*nord exactly. The number of quadrature points is calculated as:

numPoints = ceil(gaussOrder / 2) + 1

where gaussOrder = 2 * nord. For stability and predefined table limits, the number of points is saturated at 11.

Basic checks ensure that nord is non-negative and within the supported range (<=11).

Parameters:

nord – [in] Basis order driving the quadrature degree.
quadrature – [out] Rule whose numPoints is set.
nord – [in] The polynomial order of the element.
quadrature – [out] Pointer to a Quadrature1D structure where numPoints will be set.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS always.

PetscErrorCode computeNum2DQuadraturePoints(const PetscInt nord, Quadrature2D *quadrature)

Determines the number of 2D quadrature points for a triangle.

Uses a precomputed lookup table for triangular quadrature rules to determine the minimum number of points needed to exactly integrate polynomials of degree up to nord over a triangle. Valid orders are 0 through 19. Each table entry corresponds to the recommended number of quadrature points for that order. Basic checks ensure nord is within the supported range.

Parameters:

nord – [in] Basis order driving the quadrature degree.
quadrature – [out] Rule whose numPoints is set.
nord – [in] The polynomial order of the element (0 ≤ nord ≤ 19).
quadrature – [out] Pointer to a Quadrature2D structure where numPoints will be set.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS always.

PetscErrorCode computeNum3DQuadraturePoints(const PetscInt nord, Quadrature3D *quadrature)

Determines the number of Gauss quadrature points for a tetrahedron.

Computes the required number of Gauss points using a precomputed lookup table for tetrahedral quadrature rules. Supports integration orders corresponding to 1 ≤ 2*nord ≤ 12. Basic checks ensure that nord produces a valid integration order. The table maps each Gauss order directly to the recommended number of quadrature points for exact integration over a tetrahedron.

Parameters:

nord – [in] Basis order driving the quadrature degree.
quadrature – [out] Rule whose numPoints is set.
nord – [in] The basis polynomial order of the element. The required integration order is 2*nord.
quadrature – [out] Pointer to a Quadrature3D structure where numPoints will be set.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success.

PetscErrorCode compute1DQuadraturePoints(Quadrature1D *quadrature)

Populates 1D Gauss-Legendre quadrature points and weights.

The function selects the correct precomputed 1D Gauss-Legendre table based on numPoints and copies the values into the quadrature structure. Supports up to 11 points, corresponding to polynomial integration orders up to 21 (exact for polynomials of degree 2*numPoints-1).

Parameters:

quadrature – [inout] Rule (numPoints set) whose points/weights fill.
quadrature – [inout] Pointer to a Quadrature1D structure. The numPoints field must be set before calling. After execution, points and weights are filled with the appropriate Gauss-Legendre points and weights.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success.

PetscErrorCode compute2DQuadraturePoints(Quadrature2D *quadrature)

Populates 2D Gauss quadrature points and weights for a triangle.

Selects the precomputed 2D triangular quadrature table corresponding to numPoints and normalizes the weights. Supports standard Gauss rules with 1, 3, 4, 6, 7, 12, 13, 16, 19, 25, 27, 33, 37, 42, 48, 52, 61, 70, and 73 points. Exact for polynomials of degree up to the specified order.

Parameters:

quadrature – [inout] Rule (numPoints set) whose points/weights fill.
quadrature – [inout] Pointer to a Quadrature2D structure. The numPoints field must be set before calling. After execution, points and weights are filled with the appropriate 2D Gauss points and weights.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success.

PetscErrorCode compute3DQuadraturePoints(Quadrature3D *quadrature)

Populates 3D Gauss quadrature points and weights for a tetrahedron.

Selects the appropriate precomputed 3D tetrahedral quadrature table based on numPoints and calls renormalization3DGaussPoints to normalize and copy the values into the quadrature structure. Supports the following numbers of points: 1, 4, 5, 11, 14, 24, 31, 43, 53, 126, 210, which correspond to polynomial integration orders from 1 up to 12.

Parameters:

quadrature – [inout] Rule (numPoints set) whose points/weights fill.
quadrature – [inout] Pointer to a Quadrature3D structure. The numPoints field must be set before calling. After execution, points and weights are filled with the appropriate 3D Gauss points and weights.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success.

PetscErrorCode tetrahedronXYZToReference(const PetscReal coordinates[NUM_VERTICES_PER_CELL * NUM_DIMENSIONS], const PetscReal point[NUM_DIMENSIONS], PetscReal XiEtaZeta[NUM_DIMENSIONS])

Maps global Cartesian coordinates to reference-tetrahedron coordinates.

This function computes the reference coordinates of a point with respect to a physical tetrahedron. The reference tetrahedron is assumed to have vertices at (0,0,0), (1,0,0), (0,1,0), and (0,0,1).

The mapping is affine, so the inverse is computed using Cramer’s rule via triple products (determinants of 3x3 matrices). Specifically, if v0, v1, v2 are the vectors from the first vertex to the other three vertices and vp is the vector from the first vertex to the point, the reference coordinates are given by: xi = det(vp, v1, v2) / det(v0, v1, v2) eta = det(v0, vp, v2) / det(v0, v1, v2) zeta= det(v0, v1, vp) / det(v0, v1, v2)

The function checks for degenerate tetrahedra (zero or near-zero volume) and raises an error if detected.

Parameters:

coordinates – [in] Cell vertex coordinates (4 vertices × 3 dims).
point – [in] Global point to map.
XiEtaZeta – [out] Reference coordinates (ξ, η, ζ) of the point.
coordinates – [in] Array of size 12 containing the spatial coordinates of the tetrahedron’s 4 vertices in the order [x0, y0, z0, x1, y1, z1, x2, y2, z2, x3, y3, z3].
point – [in] The global Cartesian coordinates [x, y, z] of the point to transform.
XiEtaZeta – [out] Output array [xi, eta, zeta] representing the coordinates of the point in the reference tetrahedron.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS always.

PetscErrorCode computeVectorRotation(const PetscReal azimuth, const PetscReal dip, PetscReal rotationVector[NUM_DIMENSIONS])

Computes a 3D unit vector from sequential azimuth/dip rotations.

The function starts with the base vector [1, 0, 0] (aligned with the x-axis). It converts the input angles from degrees to radians. Rotations are applied sequentially:

x-y plane rotation by azimuth (matrix M1)
x-z plane rotation by dip (matrix M2)
y-z plane rotation (matrix M3), currently hardcoded with angle tetha = 0.

The overall rotation matrix is M = M1 * M2 * M3. The final rotated vector is obtained by multiplying the base vector by this rotation matrix and stored in rotationVector.

Note

The y-z plane rotation is currently fixed at 0 degrees, but the code structure allows for future generalization.

Parameters:

azimuth – [in] Azimuth angle (degrees).
dip – [in] Dip angle (degrees).
rotationVector – [out] Resulting unit direction vector.
azimuth – [in] Rotation angle in the x-y plane (degrees).
dip – [in] Rotation angle in the x-z plane (degrees).
rotationVector – [out] Output 3D vector after rotation.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS always.

PetscErrorCode computeElementalMatrices(const FEMSpace *fem, const Cell *cell, const Quadrature3D *quadrature, PetscReal **Me, PetscReal **Ke)

Computes the elemental mass and stiffness matrices for a cell.

This function performs the following steps to assemble the local elemental matrices for a tetrahedral element in 3D:

Tensor Setup
- Define the permittivity tensor e_r from the cell conductivity.
- Define the magnetic permeability tensor mu_r (currently identity).
Memory Allocation
- Allocate arrays for basis functions (Ni), their curls (NiCurl), derivatives (Dx_Ni, Dy_Ni, Dz_Ni), and Nédélec coefficients.
Reset Elemental Matrices
- Initialize Me and Ke to zero.
Nédélec Coefficients
- Compute the cell-local Nédélec basis coefficients (and, for nord=1, the monomial derivatives) via the per-order dispatch table fem->ops->computeCoefficients.
Loop Over Quadrature Points
- For each Gauss point:
  - Transform coordinates to the physical element.
  - Compute basis functions at the Gauss point via fem->ops->computeBasis.
  - Compute the mass matrix contribution: Me_ij += w_q * (N_i · (ε_r * N_j)) * sign_i * sign_j * det(J)
  - Compute the curls of the basis functions via fem->ops->computeCurls.
  - Compute the stiffness matrix contribution: Ke_ij += w_q * (curl(N_i) · μ_r * curl(N_j)) * sign_i * sign_j * det(J)
  - Edge orientation signs from cell->orientation are applied to ensure global consistency.
Cleanup
- Free all dynamically allocated arrays for Ni, NiCurl, derivatives, and coefficients.

Note

The stiffness matrix assumes isotropic magnetic permeability (μ_r = 1.0).

Parameters:

fem – [in] Finite-element space descriptor (order, DOF counts, ops).
cell – [in] Cell geometry and orientation.
quadrature – [in] 3D quadrature rule.
Me – [out] Elemental mass matrix (numDofInCell²).
Ke – [out] Elemental stiffness matrix (numDofInCell²).
nord – [in] Polynomial order of the Nédélec basis functions (1..6, dispatched via fem->ops).
numDofInCell – [in] Number of degrees of freedom per cell (edges).
cell – [in] Pointer to the Cell structure containing:
- Jacobian matrix and determinant
- Orientation of edges
- Conductivity tensor for material properties
quadrature – [in] Pointer to the Quadrature3D structure containing:
- Gauss points
- Weights for tetrahedral integration
Me – [out] Elemental mass matrix (numDofInCell x numDofInCell), representing \int_T (ε_r * N_i) · N_j dV.
Ke – [out] Elemental stiffness matrix (numDofInCell x numDofInCell), representing \int_T (μ_r^{-1} * curl(N_i)) · curl(N_j) dV.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success.

PetscErrorCode buildDofSigns(const Cell *cell, const FEMSpace *fem, PetscInt signs[])

Builds the per-DOF sign vector for a cell at the chosen Nédélec order.

nord=1: one DOF per edge; signs[e] = cell->orientation.edgeSigns[e]. nord=2: faces-first with 2 DOFs per entity. Face DOFs are +1 (canonical- geometry q-vectors agree across adjacent cells); edge-DOF signs come from cell->orientation.edgeSigns[0..5].

Builds the per-DOF sign vector for a cell at the chosen Nédélec order.

For the unified hierarchical Nédélec basis (nord = 1..6) all per-DOF sign multipliers are +1: orientation is already encoded inside the reference shape functions (OrientE / OrientTri inside shape3DETet), so multiplying by cell->orientation.edgeSigns[] here would double- apply orientation and break tangential continuity across shared faces.

Parameters:

cell – [in] Cell with computed orientation.
fem – [in] Finite-element space descriptor.
signs – [out] Caller-provided array of length fem->numDofInCell.
cell – [in] Cell (unused under the hierarchical basis; kept for API symmetry with hypothetical non-hierarchical paths).
fem – [in] FEM space - only fem->numDofInCell is consulted.
signs – [out] Array of length fem->numDofInCell; all entries set to +1 on return.

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

Returns:

PetscErrorCode PETSC_SUCCESS on success.

PetscErrorCode evaluateNedelecBasis(const FEMSpace *fem, const Cell *cell, const PetscReal point[NUM_DIMENSIONS], PetscReal **coeffs, PetscReal **Dx_Ni, PetscReal **Dy_Ni, PetscReal **Dz_Ni, PetscReal **Ni, PetscReal **NiCurl)

Evaluates the Nédélec basis (and optionally curls) at a reference point.

Dispatches on fem->nord. Pass NiCurl = NULL to skip curl evaluation.

Parameters:

fem – [in] Finite-element space descriptor.
cell – [in] Cell geometry and orientation.
point – [in] Reference-cell evaluation point.
coeffs – [inout] Basis-coefficient workspace (order-dependent).
Dx_Ni – [inout] Workspace for ∂/∂x terms (used at nord=1).
Dy_Ni – [inout] Workspace for ∂/∂y terms (used at nord=1).
Dz_Ni – [inout] Workspace for ∂/∂z terms (used at nord=1).
Ni – [out] Basis-function values at the point.
NiCurl – [out] Basis-curl values (NULL to skip).

Returns:

PetscErrorCode PETSC_SUCCESS on success, or a PetscError code otherwise.

const NedelecOps *nedelecOpsForOrder(PetscInt nord)

Returns the per-order Nédélec ops table.

Parameters:

nord – [in] Basis order.

Returns:

Pointer to the ops table, or NULL if nord is unsupported.

struct Quadrature1D

#include <hvfem.h>

1D Gauss quadrature rule (points and weights).

Public Members

PetscInt numPoints: Number of quadrature points.

PetscReal *points: Quadrature point coordinates.

PetscReal *weights: Quadrature weights.

struct Quadrature2D

#include <hvfem.h>

2D (triangle) Gauss quadrature rule (points and weights).

Public Members

PetscInt numPoints: Number of quadrature points.

PetscReal **points: Quadrature point coordinates (per point).

PetscReal *weights: Quadrature weights.

struct Quadrature3D

#include <hvfem.h>

3D (tetrahedron) Gauss quadrature rule (points and weights).

Public Members

PetscInt numPoints: Number of quadrature points.

PetscReal **points: Quadrature point coordinates (per point).

PetscReal *weights: Quadrature weights.

struct NedelecOps

#include <hvfem.h>

Per-order Nédélec dispatch table.

Each supported basis order registers one of these; hot paths call through the table instead of a switch (nord). Adding a new order means: define the four per-order static helpers in hvfem.c, build a NedelecOps instance for them, and register it in nedelecOpsForOrder().

Uniform signatures - some fields are unused for a given order:

computeCoefficients: nord=1 fills coeffs AND Dx/Dy/Dz; nord>=2 only fills coeffs (Dx/Dy/Dz still passed but untouched).
computeBasis : evaluates Ni at a reference-cell point.
computeCurls : evaluates NiCurl at a reference-cell point. nord=1 uses Dx/Dy/Dz; nord>=2 uses coeffs+point.
buildGradientMatrix: TOPOLOGICAL discrete gradient G against P1 H1 (numDofInCell x NUM_H1_DOF_PER_CELL), used as the BDDC Nédélec hint via PCBDDCSetDiscreteGradient. Only the lowest-order Whitney row per mesh edge is nonzero (±1 vertex incidence); K·G is NOT zero at nord >= 2.

Public Members

PetscErrorCode (*computeCoefficients)(const Cell *cell, PetscReal **coeffs, PetscReal **Dx_Ni, PetscReal **Dy_Ni, PetscReal **Dz_Ni): Computes basis coefficients (and, at nord=1, derivative tables).

PetscErrorCode (*computeBasis)(const Cell *cell, const PetscReal point[NUM_DIMENSIONS], const PetscReal *const *coeffs, PetscReal **Ni): Evaluates basis values Ni at a reference-cell point.

PetscErrorCode (*computeCurls)(const Cell *cell, const PetscReal point[NUM_DIMENSIONS], const PetscReal *const *coeffs, const PetscReal *const *Dx_Ni, const PetscReal *const *Dy_Ni, const PetscReal *const *Dz_Ni, PetscReal **NiCurl): Evaluates basis curls NiCurl at a reference-cell point.

PetscErrorCode (*buildGradientMatrix)(const FEMSpace *fem, const Cell *cell, const Quadrature1D *quadrature1d, PetscReal **gradientMatrix): Builds the topological discrete-gradient matrix (BDDC hint).