Directional C++ Library

Eigen is a C++ template library designed for linear algebra operations. It provides high-performance implementations for:

• Matrices and Vectors: Core data structures for linear algebra.
• Numerical Solvers: Algorithms for solving linear systems and other numerical problems.
• Related Algorithms: A wide range of mathematical algorithms built on top of the matrix/vector primitives.

Because it is a template library, it is designed for high performance and can be integrated into C++ projects to handle complex mathematical computations.

Most tensor operations return an <Operation> (a deferred expression) rather than a result. You cannot call .dimensions() or operator() directly on these expressions.

To work with a deferred operation as if it were a real tensor, wrap it in a TensorRef. This is particularly useful in templated code where you want to handle both raw Tensor objects and slices/expressions agnostically.

Workflow:

1. Perform an operation (returns <Operation>).
2. Wrap the result in a TensorRef.
3. Use TensorRef to access dimensions or evaluate the expression into a concrete tensor.

// If 'op' is a deferred operation (e.g. from a slice or math op):
// Eigen::TensorRef<Eigen::Tensor<float, 2>> ref = op;
// Now you can use ref.dimensions() or ref(i, j)

You can reduce the degrees of freedom in the integrated functions by enforcing linear relations $F = U imes f$. This is done by providing a matrix $U$ to the linRed field in IntegrationData.

Commonly used reductions include:

• Triangular symmetry: Use set_triangular_symmetry(int N).
• Sign symmetry: Automatically assumed when $N$ is even ($U = [Id_{n imes n}; -Id_{n imes n}]$).

Warning: The reduction must commute with the permutation around singularities to avoid algorithm failure.

Hodge decomposition splits an $i$-form $z_i$ into three unique components: an exact component ($d_{i-1}f_{i-1}$), a coexact component ($\delta_i g_{i+1}$), and a harmonic component ($h_i$).

To perform this decomposition, use the following logic (conceptually implemented via Hodge_decomposition()):

1. Exact component: Use project_exact(d, M, cochain, prevCochain, exactCochain). If the operator $d$ has a trivial null-space, it solves the Poisson equation. If it has a non-trivial null-space, use the GagueFixing flag to solve for the minimum $2$-norm solution.
2. Coexact component: Use project_coexact(dNext, M, MNext, cochain, k, coexactCochain, nextCochain).
3. Harmonic component: Calculated as harmCochain = cochain - exactCochain - coexactCochain.

All Tensor methods (except for direct assignments or specific evaluation calls) return non-evaluated Operations. These operations can be chained together to form complex expressions.

Crucially, the chain of operations is evaluated lazily. The actual computation typically only occurs when the operation is assigned to a tensor. This allows the library to optimize the expression before execution.

Directional provides flexibility in discretization by allowing you to compute fields on either a vertex-based or a face-based tangent bundle.

Available types:

• directional::IntrinsicVertexTangentBundle vtb (Vertex-based)
• directional::PCFaceTangentBundle ftb (Face-based)

To visualize both, initialize both bundles with the mesh, set them as Cartesian fields with unique names and ordinal numbers, and use viewer.toggle_cartesian_field(ordinal, name) or viewer.toggle_singularities(ordinal, name) to switch between them in the UI.

directional::IntrinsicVertexTangentBundle vtb;
directional::PCFaceTangentBundle ftb;

// ... setup mesh ...

ftb.init(mesh);
vtb.init(mesh);

// ... set fields with unique ordinals ...
viewer.set_cartesian_field(rawFaceField, "Face-Based Field", 0);
viewer.set_cartesian_field(rawVertexField, "Vertex-Based Field", 1);
viewer.launch();

When working with Eigen Tensors, the documentation uses several pseudo-types to describe methods and operations:

• <Tensor-Type>::Dimensions: An array-like object of ints representing the tensor's dimensions. It has a size attribute (the rank) and can be indexed like an array.
• <Tensor-Type>::Index: An integer type used for indexing along dimensions (returned by dimension(), size(), or used in operator()).
• <Tensor-Type>::Scalar: The underlying data type of the tensor elements (e.g., float for Tensor<float, 2>).
• <Operation>: Indicates that a method returns a deferred operation rather than a concrete tensor. An <Operation> must be evaluated (e.g., by assigning it to a tensor) before its values can be accessed directly.

Eigen Tensor operations (like +, *, .exp()) do not perform computations immediately. Instead, they construct a lightweight "tensor operator" object (e.g., TensorCwiseBinaryOp). This allows for lazy evaluation, where a complex tree of operations is built and only executed when the result is assigned to a concrete Tensor type. This mechanism enables significant compiler optimizations and performance gains.

Example of building an expression tree:

// This constructs a tree of operators, it does NOT perform addition.
auto t4 = t1 + t2;

// This constructs a nested tree: TensorCwiseBinaryOp(t1, TensorCwiseUnaryOp(t2, 0.3f))
auto t5 = t1 + t2 * 0.3f;

Tensor<float, 3> t1(2, 3, 4);
Tensor<float, 3> t2(2, 3, 4);
// Set t3 to the element wise sum of t1 and t2
Tensor<float, 3> t3 = t1 + t2;

In the Eigen Tensor library, scalar values are often represented as tensors with a rank of 0 and a size of 1.

For example:

• Tensor<T, N>::maximum() returns a Tensor<T, 0>.
• The inner product of two 1D tensors (via contractions) returns a 0D tensor.