Algorithmic Functions

Overview

The algorithmic module provides advanced numerical algorithms that operate on meshes. The primary algorithm is HarmonicDiffusion, which solves the Laplace equation on triangulated surfaces using iterative relaxation. This is useful for smoothly interpolating sparse constraints across a mesh, such as boundary conditions, temperature fields, or displacement constraints.

HarmonicDiffusion

HarmonicDiffusion Class


template <size_t N>
class HarmonicDiffusion {
public:
    // Constructor: takes a mesh and an initial value for all nodes
    HarmonicDiffusion(const Mesh<N> &mesh, double init_value = 0.0);

    // Solver parameters
    void setMaxIter(size_t max_iterations);   // Default: 10000
    void setEps(double epsilon);              // Convergence tolerance (default: 1e-6)
    void setEpsilon(double epsilon);          // Alias for setEps

    // Add a constraint: fix the value at a specific vertex
    void addConstraint(uint32_t vertex_index, double value);

    // Constrain all border nodes to a given value
    void constrainBorders(double value);

    // Solve the Laplace equation and return the solution
    Serie<double> solve();
};

The harmonic diffusion solver iteratively relaxes the solution toward the Laplace equation (zero second derivative) on the mesh surface. At each iteration, each unconstrained node is updated to the average of its neighbors' values. The process converges when the maximum change between iterations falls below the specified tolerance (epsilon).

Basic Usage

Harmonic Diffusion Example


#include <dataframe/geo/mesh.h>
#include <dataframe/algo/harmonic_diffusion.h>

// Create or load a mesh
df::Mesh3D mesh(vertices, triangles);

// Create solver with initial value 0
df::HarmonicDiffusion<3> solver(mesh, 0.0);

// Set solver parameters
solver.setMaxIter(5000);
solver.setEps(1e-8);

// Fix values at specific vertices
solver.addConstraint(0, 100.0);    // Vertex 0 = 100
solver.addConstraint(50, 0.0);     // Vertex 50 = 0
solver.addConstraint(100, 50.0);   // Vertex 100 = 50

// Solve
auto result = solver.solve();
// result is a Serie<double> with smoothly interpolated values

Border Constraints

Constraining Borders Example


// Heat diffusion: fix borders to 0, internal source to 100
df::HarmonicDiffusion<3> heat_solver(mesh, 0.0);

// Fix all border nodes to 0 (cold boundary)
heat_solver.constrainBorders(0.0);

// Add a hot spot in the interior
heat_solver.addConstraint(center_vertex_id, 100.0);

heat_solver.setMaxIter(10000);
heat_solver.setEps(1e-6);

auto temperature = heat_solver.solve();

// Add result as mesh attribute
mesh.addVertexAttribute("temperature", temperature);

Complete Example

Parameterization on a Surface


#include <dataframe/Serie.h>
#include <dataframe/geo/mesh.h>
#include <dataframe/geo/geometry.h>
#include <dataframe/algo/harmonic_diffusion.h>
#include <iostream>

int main() {
    // Generate a sphere mesh
    auto [verts, tris] = df::geo::generateSphere(4, 1.0);
    df::Mesh3D sphere(verts, tris);

    // Create harmonic diffusion solver
    df::HarmonicDiffusion<3> solver(sphere, 0.5);
    solver.setMaxIter(20000);
    solver.setEps(1e-8);

    // Find the north and south pole vertices (max/min z)
    size_t north_pole = 0, south_pole = 0;
    double max_z = -1e30, min_z = 1e30;
    for (size_t i = 0; i < verts.size(); ++i) {
        if (verts[i][2] > max_z) { max_z = verts[i][2]; north_pole = i; }
        if (verts[i][2] < min_z) { min_z = verts[i][2]; south_pole = i; }
    }

    // Constrain poles
    solver.addConstraint(north_pole, 1.0);
    solver.addConstraint(south_pole, 0.0);

    // Solve for a smooth scalar field
    auto field = solver.solve();
    sphere.addVertexAttribute("latitude_param", field);

    std::cout << "Solved harmonic diffusion on sphere with "
              << sphere.vertexCount() << " vertices\n";
    std::cout << "North pole value: " << field[north_pole] << "\n";
    std::cout << "South pole value: " << field[south_pole] << "\n";

    return 0;
}

Implementation Notes

The solver uses Jacobi iteration on the graph Laplacian of the mesh.
Convergence speed depends on mesh connectivity and the distribution of constraints.
For large meshes, consider increasing maxIter and monitoring convergence.
Constrained vertices are never modified during iteration.
Works with both 2D and 3D meshes (Mesh<2> and Mesh<3>).