cellsurfaceinterp.h

Go to the documentation of this file.

 
#ifndef CELL_SURFACE_INTERP_H
#define CELL_SURFACE_INTERP_H
 
#include <QPointF>
 
#include <algorithm>
 
namespace CellSurfaceInterp
{
 
inline double depthAt(const QPointF &p,
                      const QPointF &a, const QPointF &b, const QPointF &c,
                      double z0, double z1, double z2,
                      double sd0, double sd1, double sd2,
                      bool *degenerate = nullptr)
{
    if (degenerate) *degenerate = false;
 
    const double v0x = c.x() - a.x(), v0y = c.y() - a.y();
    const double v1x = b.x() - a.x(), v1y = b.y() - a.y();
    const double v2x = p.x() - a.x(), v2y = p.y() - a.y();
    const double d00 = v0x * v0x + v0y * v0y;
    const double d01 = v0x * v1x + v0y * v1y;
    const double d11 = v1x * v1x + v1y * v1y;
    const double d20 = v2x * v0x + v2y * v0y;
    const double d21 = v2x * v1x + v2y * v1y;
    const double denom = d00 * d11 - d01 * d01;
    if (denom == 0.0) {
        if (degenerate) *degenerate = true;
        return 0.0;
    }
    const double u = (d11 * d20 - d01 * d21) / denom;   // weight for c
    const double v = (d00 * d21 - d01 * d20) / denom;   // weight for b
    const double w = 1.0 - u - v;                       // weight for a
 
    // Only WET corners (sd > 0 — water actually stands above their bed)
    // SUPPLY surface. The other two states carry no water at their corner:
    //   sd == 0  no data — η undefined;
    //   sd  < 0  η valid but below the corner's own bed — informational only.
    // A below-bed corner must neither drag the blend down (a thin flank film
    // pooled at a wall base would notch the pool surface toward the wall)
    // nor raise the driving-head cap (a high bank corner with η above the
    // pool would license the blend to climb the wall face) — both were the
    // residual profile artifacts after the first iteration of this fix.
    const bool wet[3] = { sd0 > 0.0, sd1 > 0.0, sd2 > 0.0 };
    const int nWet = int(wet[0]) + int(wet[1]) + int(wet[2]);
    // Vertex-scoped dryness: no wet corner → nothing can flood this cell.
    if (nWet == 0) return 0.0;
 
    const double eta[3] = { z0 + sd0, z1 + sd1, z2 + sd2 };
    const double groundInterp = w * z0 + v * z1 + u * z2;
 
    // Driving-head cap over the WET corners' η — the head the water actually
    // supplies; guards extrapolating weights at/outside the edge.
    bool   any    = false;
    double maxEta = 0.0;
    for (int k = 0; k < 3; ++k) {
        if (!wet[k]) continue;
        if (!any || eta[k] > maxEta) { maxEta = eta[k]; any = true; }
    }
 
    double blend;
    if (nWet == 3) {
        // All corners wet — the signed-depth blend equals η_blend − ground
        // algebraically; computed in sd space so the fully-wet case stays
        // bit-identical to the pre-extrapolation implementation.
        blend = w * sd0 + v * sd1 + u * sd2;
    } else {
        // Renormalized wet-corner blend: η is the weight-renormalized average
        // of the wet corners' η. On an all-wet edge this equals the plain
        // barycentric lerp (continuous with fully-wet neighbours); toward a
        // non-supplying corner the weights shrink together, so η holds level
        // (zero gradient in the dry direction) and depth = max(0, η − ground)
        // lands the waterline at the sub-cell bed intercept — the surface
        // extends flat to the wall face instead of stopping at a cell edge.
        double wetW = 0.0, wetE = 0.0;
        const double wgt[3] = { w, v, u };
        for (int k = 0; k < 3; ++k) {
            if (!wet[k]) continue;
            wetW += wgt[k];
            wetE += wgt[k] * eta[k];
        }
        // Degenerate weight mass (sample at/beyond the dry corner, where the
        // wet weights cancel): fall back to the driving head — the cap and
        // the 0-floor still confine the result to the physical band.
        const double etaS = (wetW > 1e-12) ? wetE / wetW : maxEta;
        blend = etaS - groundInterp;
    }
 
    const double capDepth = maxEta - groundInterp;
    return std::max(0.0, std::min(blend, capDepth));
}
 
} // namespace CellSurfaceInterp
 
#endif // CELL_SURFACE_INTERP_H