ImplicitSurfaceNetwork/surface_integral/interface/mesh_algorithm.hpp

#include <cmath>
#include <cstdio>
#include <solve.h>  // polymesh_t


// Vector subtraction: a - b
inline vector3d vec3_sub(const vector3d& a, const vector3d& b) {
    return { a.x - b.x, a.y - b.y, a.z - b.z };
}

// Vector cross product: a × b
inline vector3d vec3_cross(const vector3d& a, const vector3d& b) {
    return {
        a.y * b.z - a.z * b.y,
        a.z * b.x - a.x * b.z,
        a.x * b.y - a.y * b.x
    };
}

// Vector dot product: a · b
inline double vec3_dot(const vector3d& a, const vector3d& b) {
    return a.x * b.x + a.y * b.y + a.z * b.z;
}

// Vector length: ||a||
inline double vec3_norm(const vector3d& a) {
    return std::sqrt(a.x*a.x + a.y*a.y + a.z*a.z);
}

// Compute triangle area: 0.5 * || (b-a) × (c-a) ||
inline double triangle_area(const vector3d& a, const vector3d& b, const vector3d& c) {
    vector3d ab = vec3_sub(b, a);
    vector3d ac = vec3_sub(c, a);
    vector3d cross_product = vec3_cross(ab, ac);
    return 0.5 * vec3_norm(cross_product);
}

// Compute triangle signed volume contribution (w.r.t origin): (a · (b × c)) / 6
inline double triangle_signed_volume(const vector3d& a, const vector3d& b, const vector3d& c) {
    vector3d cross_product = vec3_cross(b, c);
    return vec3_dot(a, cross_product) / 6.0;
}

// Compute polygon face area (triangulate as fan)
double polygon_area(const vector3d* vertices, const uint32_t* face_indices, uint32_t n) {
    if (n < 3) return 0.0;
    double area = 0.0;
    const vector3d& v0 = vertices[face_indices[0]];
    for (uint32_t i = 1; i < n - 1; ++i) {
        const vector3d& v1 = vertices[face_indices[i]];
        const vector3d& v2 = vertices[face_indices[i + 1]];
        area += triangle_area(v0, v1, v2);
    }
    return area;
}

// Compute polygon signed volume contribution (sum of triangle contributions)
double polygon_signed_volume(const vector3d* vertices, const uint32_t* face_indices, uint32_t n) {
    if (n < 3) return 0.0;
    double volume = 0.0;
    const vector3d& v0 = vertices[face_indices[0]];
    for (uint32_t i = 1; i < n - 1; ++i) {
        const vector3d& v1 = vertices[face_indices[i]];
        const vector3d& v2 = vertices[face_indices[i + 1]];
        volume += triangle_signed_volume(v0, v1, v2);
    }
    return volume;
}

// 🟩 Main: compute total mesh surface area
double compute_surface_area(const polymesh_t* mesh) {
    if (!mesh || !mesh->vertices || !mesh->faces || !mesh->vertex_counts) {
        return 0.0;
    }

    double total_area = 0.0;
    const uint32_t* face_ptr = mesh->faces;  // current index into faces array

    for (uint32_t i = 0; i < mesh->num_faces; ++i) {
        uint32_t n = mesh->vertex_counts[i];
        total_area += polygon_area(mesh->vertices, face_ptr, n);
        face_ptr += n;  // move to the next face's start index
    }

    return total_area;
}

// 🟦 Main: compute total mesh volume (requires closed mesh with outward normals)
double compute_volume(const polymesh_t* mesh) {
    if (!mesh || !mesh->vertices || !mesh->faces || !mesh->vertex_counts) {
        return 0.0;
    }

    double total_volume = 0.0;
    const uint32_t* face_ptr = mesh->faces;

    for (uint32_t i = 0; i < mesh->num_faces; ++i) {
        uint32_t n = mesh->vertex_counts[i];
        total_volume += polygon_signed_volume(mesh->vertices, face_ptr, n);
        face_ptr += n;
    }

    return std::abs(total_volume);  // ensure positive volume
}