meshcellgeom.h

Go to the documentation of this file.

 
#ifndef OPENSWMMVIS_MESH_MESHCELLGEOM_H
#define OPENSWMMVIS_MESH_MESHCELLGEOM_H
 
#include "mesh/meshresult.h"
 
#include <QPointF>
#include <QVector>
 
#include <array>
#include <cmath>
 
namespace mesh {
 
constexpr int kEdgeStride = 4;
 
constexpr int edgeSlot(int cell, int e) noexcept { return cell * kEdgeStride + e; }
constexpr int slotCell(int slot) noexcept { return slot / kEdgeStride; }
constexpr int slotLocal(int slot) noexcept { return slot % kEdgeStride; }
constexpr int edgeSlotCount(int nCells) noexcept { return nCells * kEdgeStride; }
 
inline void edgeEndpoints(const MeshTriangle &t, int k, int &a, int &b) noexcept
{
    const int nv = t.vertexCount();
    a = t.vertex((k + 1) % nv);
    b = t.vertex((k + 2) % nv);
}
 
inline bool rt0CellDischarge(int cell, int nv,
                             const float *flux, const float *length,
                             const float *nx, const float *ny,
                             double &qx, double &qy) noexcept
{
    constexpr double kQMax = 10.0;
    double a00 = 0.0, a01 = 0.0, a11 = 0.0;  // NᵀN entries
    double b0  = 0.0, b1  = 0.0;             // Nᵀ b entries
    for (int e = 0; e < nv; ++e) {
        const int idx = edgeSlot(cell, e);
        const double len = length[idx];
        if (len <= 1e-12) continue;
        double q = double(flux[idx]) / len;
        if (!std::isfinite(q)) continue;
        if (q >  kQMax) q =  kQMax;
        if (q < -kQMax) q = -kQMax;
        const double enx = nx[idx], eny = ny[idx];
        a00 += enx * enx;
        a01 += enx * eny;
        a11 += eny * eny;
        b0  += enx * q;
        b1  += eny * q;
    }
    const double det = a00 * a11 - a01 * a01;
    // NaN-robust degeneracy gate: `abs(NaN) < eps` is false, so the inverted
    // form is required to reject the cell instead of emitting NaN.
    if (!(std::abs(det) >= 1e-12)) return false;
    const double inv_det = 1.0 / det;
    qx = ( a11 * b0 - a01 * b1) * inv_det;
    qy = (-a01 * b0 + a00 * b1) * inv_det;
    return true;
}
 
struct CellGeom
{
    double area = 0.0;                 
    QPointF centroid;                  
    double zMean = 0.0;                
    int    nSub = 1;                   
    std::array<std::array<int, 3>, 2> sub{};   
    int    diagCase = 0;               
};
 
inline CellGeom cellGeom(const QVector<MeshVertex> &vertices, const MeshTriangle &t)
{
    CellGeom g;
    auto triArea = [&](int a, int b, int c) {
        const QPointF &A = vertices[a].xy, &B = vertices[b].xy, &C = vertices[c].xy;
        return 0.5 * std::abs((B.x() - A.x()) * (C.y() - A.y()) - (C.x() - A.x()) * (B.y() - A.y()));
    };
    auto triCentroid = [&](int a, int b, int c) {
        const QPointF &A = vertices[a].xy, &B = vertices[b].xy, &C = vertices[c].xy;
        return QPointF((A.x() + B.x() + C.x()) / 3.0, (A.y() + B.y() + C.y()) / 3.0);
    };
    if (!t.isQuad()) {
        g.nSub = 1;
        g.sub[0] = {t.v0, t.v1, t.v2};
        g.area = triArea(t.v0, t.v1, t.v2);
        g.centroid = triCentroid(t.v0, t.v1, t.v2);
        g.zMean = (vertices[t.v0].z + vertices[t.v1].z + vertices[t.v2].z) / 3.0;
        return g;
    }
    const int v[4] = {t.v0, t.v1, t.v2, t.v3};
    // Elevation ordering of the cyclic positions (stable insertion sort).
    int p[4] = {0, 1, 2, 3};
    for (int i = 1; i < 4; ++i) {
        const int key = p[i];
        int j = i - 1;
        while (j >= 0 && vertices[v[p[j]]].z > vertices[v[key]].z) { p[j + 1] = p[j]; --j; }
        p[j + 1] = key;
    }
    const int n1 = p[0], n2 = p[1], n3 = p[2], n4 = p[3];
    auto adjacent = [](int a, int b) { const int d = (a - b + 4) % 4; return d == 1 || d == 3; };
    int t1[3], t2[3];
    if (!adjacent(n1, n4))      { g.diagCase = 1; t1[0]=n1; t1[1]=n2; t1[2]=n4; t2[0]=n1; t2[1]=n3; t2[2]=n4; }
    else if (adjacent(n2, n1))  { g.diagCase = 2; t1[0]=n1; t1[1]=n2; t1[2]=n4; t2[0]=n2; t2[1]=n3; t2[2]=n4; }
    else                        { g.diagCase = 3; t1[0]=n1; t1[1]=n3; t1[2]=n4; t2[0]=n2; t2[1]=n3; t2[2]=n4; }
    g.nSub = 2;
    g.sub[0] = {v[t1[0]], v[t1[1]], v[t1[2]]};
    g.sub[1] = {v[t2[0]], v[t2[1]], v[t2[2]]};
    const double a1 = triArea(g.sub[0][0], g.sub[0][1], g.sub[0][2]);
    const double a2 = triArea(g.sub[1][0], g.sub[1][1], g.sub[1][2]);
    g.area = a1 + a2;
    const QPointF c1 = triCentroid(g.sub[0][0], g.sub[0][1], g.sub[0][2]);
    const QPointF c2 = triCentroid(g.sub[1][0], g.sub[1][1], g.sub[1][2]);
    g.centroid = (g.area > 0.0)
        ? QPointF((a1 * c1.x() + a2 * c2.x()) / g.area, (a1 * c1.y() + a2 * c2.y()) / g.area)
        : QPointF((vertices[t.v0].xy + vertices[t.v1].xy + vertices[t.v2].xy + vertices[t.v3].xy) / 4.0);
    g.zMean = (vertices[t.v0].z + vertices[t.v1].z + vertices[t.v2].z + vertices[t.v3].z) / 4.0;
    return g;
}
 
inline double cellSignedArea(const QVector<MeshVertex> &vertices, const MeshTriangle &t) noexcept
{
    const int nv = t.vertexCount();
    double s = 0.0;
    for (int k = 0; k < nv; ++k) {
        const QPointF &a = vertices[t.vertex(k)].xy, &b = vertices[t.vertex((k + 1) % nv)].xy;
        s += a.x() * b.y() - b.x() * a.y();
    }
    return 0.5 * s;
}
 
inline bool cellIsConvex(const QVector<MeshVertex> &vertices, const MeshTriangle &t) noexcept
{
    if (!t.isQuad()) return true;
    int sign = 0;
    for (int k = 0; k < 4; ++k) {
        const QPointF &p = vertices[t.vertex(k)].xy;
        const QPointF &q = vertices[t.vertex((k + 1) % 4)].xy;
        const QPointF &r = vertices[t.vertex((k + 2) % 4)].xy;
        const double cr = (q.x() - p.x()) * (r.y() - q.y()) - (q.y() - p.y()) * (r.x() - q.x());
        const int s = (cr > 0.0) ? 1 : (cr < 0.0) ? -1 : 0;
        if (s == 0 || (sign != 0 && s != sign)) return false;
        sign = s;
    }
    return true;
}
 
inline bool cellContains(const QVector<MeshVertex> &vertices, const MeshTriangle &t,
                         const QPointF &p, double eps = 0.0) noexcept
{
    const CellGeom g = cellGeom(vertices, t);
    for (int s = 0; s < g.nSub; ++s) {
        const QPointF &a = vertices[g.sub[s][0]].xy, &b = vertices[g.sub[s][1]].xy,
                      &c = vertices[g.sub[s][2]].xy;
        const double d = (b.y() - c.y()) * (a.x() - c.x()) + (c.x() - b.x()) * (a.y() - c.y());
        if (std::abs(d) < 1e-300) continue;
        const double l1 = ((b.y() - c.y()) * (p.x() - c.x()) + (c.x() - b.x()) * (p.y() - c.y())) / d;
        const double l2 = ((c.y() - a.y()) * (p.x() - c.x()) + (a.x() - c.x()) * (p.y() - c.y())) / d;
        const double l3 = 1.0 - l1 - l2;
        if (l1 >= -eps && l2 >= -eps && l3 >= -eps) return true;
    }
    return false;
}
 
} // namespace mesh
 
#endif // OPENSWMMVIS_MESH_MESHCELLGEOM_H