SectionGeometry.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: Apache-2.0
//
// Copyright 2026 Caleb Buahin
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

 
#ifndef OPENSWMM_ENGINE_FV_SECTION_GEOMETRY_HPP
#define OPENSWMM_ENGINE_FV_SECTION_GEOMETRY_HPP
 
#include <algorithm>
#include <cmath>
 
#include "../XSectBatch.hpp"
 
namespace openswmm::fv::secgeom {
 
inline constexpr double kPi = 3.14159265358979323846;

inline constexpr double kRectAlfMax = 0.97;
 
namespace detail {

inline double circTheta(double r, double y) {
    double c = 1.0 - y / r;
    if (c > 1.0)
        c = 1.0;
    if (c < -1.0)
        c = -1.0;
    return 2.0 * std::acos(c);
}

inline double circArea(double r, double y) {
    if (y <= 0.0)
        return 0.0;
    if (y >= 2.0 * r)
        return kPi * r * r;
    const double th = circTheta(r, y);
    return 0.5 * r * r * (th - std::sin(th));
}

inline double circWidth(double r, double y) {
    if (y <= 0.0 || y >= 2.0 * r)
        return 0.0;
    const double s = y * (2.0 * r - y);
    return (s > 0.0) ? 2.0 * std::sqrt(s) : 0.0;
}

inline double circPerim(double r, double y) {
    if (y <= 0.0)
        return 0.0;
    if (y >= 2.0 * r)
        return 2.0 * kPi * r;
    return r * circTheta(r, y);
}

template <class F>
double simpsonStep(F& f, double a, double b, double fa, double fm, double fb, double whole,
                   double eps, int depth) {
    const double m = 0.5 * (a + b);
    const double lm = 0.5 * (a + m), rm = 0.5 * (m + b);
    const double flm = f(lm), frm = f(rm);
    const double left = (m - a) / 6.0 * (fa + 4.0 * flm + fm);
    const double right = (b - m) / 6.0 * (fm + 4.0 * frm + fb);
    const double delta = left + right - whole;
    if (depth <= 0 || std::fabs(delta) <= 15.0 * eps)
        return left + right + delta / 15.0;
    return simpsonStep(f, a, m, fa, flm, fm, left, 0.5 * eps, depth - 1) +
           simpsonStep(f, m, b, fm, frm, fb, right, 0.5 * eps, depth - 1);
}

template <class F> double integrate(F f, double a, double b, double eps = 1.0e-12, int depth = 40) {
    if (!(b > a))
        return 0.0;
    const double fa = f(a), fb = f(b), fm = f(0.5 * (a + b));
    const double whole = (b - a) / 6.0 * (fa + 4.0 * fm + fb);
    return simpsonStep(f, a, b, fa, fm, fb, whole, eps, depth);
}

inline double ellipsePerim(double ax, double by, double y) {
    if (y <= 0.0)
        return 0.0;
    double c = 1.0 - y / by;
    if (c < -1.0)
        c = -1.0;
    if (c > 1.0)
        c = 1.0;
    const double phi_max = std::acos(c);
    auto ds = [&](double phi) {
        const double cx = ax * std::cos(phi), sy = by * std::sin(phi);
        return std::sqrt(cx * cx + sy * sy);
    };
    return 2.0 * integrate(ds, 0.0, phi_max, 1.0e-12 * (ax + by));
}
 
} // namespace detail

inline bool hasClosedForm(int type) {
    switch (static_cast<XSectShape>(type)) {
    case XSectShape::CIRCULAR:
    case XSectShape::FORCE_MAIN:
    case XSectShape::FILLED_CIRCULAR:
    case XSectShape::RECT_CLOSED:
    case XSectShape::RECT_OPEN:
    case XSectShape::TRAPEZOIDAL:
    case XSectShape::TRIANGULAR:
    case XSectShape::PARABOLIC:
    case XSectShape::POWERFUNC:
        return true;
    default:
        return false;
    }
}

inline bool isPolynomialOpen(int type) {
    switch (static_cast<XSectShape>(type)) {
    case XSectShape::RECT_OPEN:
    case XSectShape::TRAPEZOIDAL:
    case XSectShape::TRIANGULAR:
        return true;
    default:
        return false;
    }
}

inline double areaOfDepth(const XSectParams& xs, double h) {
    if (h <= 0.0)
        return 0.0;
    if (h > xs.y_full)
        h = xs.y_full;
    switch (static_cast<XSectShape>(xs.type)) {
    case XSectShape::CIRCULAR:
    case XSectShape::FORCE_MAIN:
        return detail::circArea(0.5 * xs.y_full, h);
    case XSectShape::FILLED_CIRCULAR: {
        const double r = 0.5 * (xs.y_full + xs.y_bot);
        return detail::circArea(r, xs.y_bot + h) - detail::circArea(r, xs.y_bot);
    }
    case XSectShape::RECT_CLOSED:
    case XSectShape::RECT_OPEN:
        return h * xs.w_max;
    case XSectShape::TRAPEZOIDAL:
        return (xs.y_bot + xs.s_bot * h) * h;
    case XSectShape::TRIANGULAR:
        return xs.s_bot * h * h;
    case XSectShape::PARABOLIC:
        return (4.0 / 3.0) * xs.r_bot * h * std::sqrt(h);
    case XSectShape::POWERFUNC:
        return xs.r_bot * std::pow(h, xs.s_bot + 1.0);
    default:
        return xsect::getAofY(xs, h);
    }
}

inline double widthOfDepth(const XSectParams& xs, double h) {
    if (h < 0.0)
        return 0.0;
    if (h > xs.y_full)
        h = xs.y_full;
    switch (static_cast<XSectShape>(xs.type)) {
    case XSectShape::CIRCULAR:
    case XSectShape::FORCE_MAIN:
        return detail::circWidth(0.5 * xs.y_full, h);
    case XSectShape::FILLED_CIRCULAR:
        return detail::circWidth(0.5 * (xs.y_full + xs.y_bot), xs.y_bot + h);
    case XSectShape::RECT_CLOSED:
    case XSectShape::RECT_OPEN:
        return xs.w_max;
    case XSectShape::TRAPEZOIDAL:
        return xs.y_bot + 2.0 * xs.s_bot * h;
    case XSectShape::TRIANGULAR:
        return 2.0 * xs.s_bot * h;
    case XSectShape::PARABOLIC:
        return 2.0 * xs.r_bot * std::sqrt(h);
    case XSectShape::POWERFUNC:
        return (h > 0.0) ? (xs.s_bot + 1.0) * xs.r_bot * std::pow(h, xs.s_bot) : 0.0;
    default:
        return xsect::getWofY(xs, h);
    }
}

inline double perimeterOfDepth(const XSectParams& xs, double h) {
    if (h <= 0.0)
        return 0.0;
    if (h > xs.y_full)
        h = xs.y_full;
    switch (static_cast<XSectShape>(xs.type)) {
    case XSectShape::CIRCULAR:
    case XSectShape::FORCE_MAIN:
        return detail::circPerim(0.5 * xs.y_full, h);
    case XSectShape::FILLED_CIRCULAR: {
        const double r = 0.5 * (xs.y_full + xs.y_bot);
        // Arc above the sediment plus the flat sediment surface.
        return detail::circPerim(r, xs.y_bot + h) - detail::circPerim(r, xs.y_bot) +
               detail::circWidth(r, xs.y_bot);
    }
    case XSectShape::RECT_CLOSED: {
        double p = xs.w_max + 2.0 * h;
        const double alf = (xs.a_full > 0.0) ? (h * xs.w_max) / xs.a_full : 0.0;
        if (alf > kRectAlfMax)
            p += (alf - kRectAlfMax) / (1.0 - kRectAlfMax) * xs.w_max;
        return p;
    }
    case XSectShape::RECT_OPEN:
        return xs.w_max + (2.0 - xs.s_bot) * h;
    case XSectShape::TRAPEZOIDAL:
        return xs.y_bot + h * xs.r_bot;
    case XSectShape::TRIANGULAR:
        return 2.0 * h * xs.r_bot;
    case XSectShape::PARABOLIC: {
        const double x = 2.0 * std::sqrt(h) / xs.r_bot;
        const double t = std::sqrt(1.0 + x * x);
        return 0.5 * xs.r_bot * xs.r_bot * (x * t + std::log(x + t));
    }
    case XSectShape::POWERFUNC: {
        // Half-width x(y) = hc·y^m; arc length 2∫√(1 + x'(y)²) dy. The
        // integrand has an integrable singularity at y → 0 when m < 1;
        // adaptive Simpson from a tiny offset handles it at build time.
        const double hc = (xs.s_bot + 1.0) * xs.r_bot / 2.0;
        const double m = xs.s_bot;
        auto ds = [&](double y) {
            const double dxdy = (y > 0.0) ? hc * m * std::pow(y, m - 1.0) : 0.0;
            return std::sqrt(1.0 + dxdy * dxdy);
        };
        const double y0 = 1.0e-9 * h;
        return 2.0 *
               (detail::integrate(ds, y0, h, 1.0e-12 * h) + std::hypot(y0, hc * std::pow(y0, m)));
    }
    default: {
        const double r = xsect::getRofY(xs, h);
        return (r > 0.0) ? xsect::getAofY(xs, h) / r : 0.0;
    }
    }
}

inline double hydRadOfDepth(const XSectParams& xs, double h) {
    if (h <= 0.0)
        return 0.0;
    if (!hasClosedForm(xs.type))
        return xsect::getRofY(xs, std::min(h, xs.y_full));
    const double p = perimeterOfDepth(xs, h);
    return (p > 0.0) ? areaOfDepth(xs, h) / p : 0.0;
}
 
} // namespace openswmm::fv::secgeom
 
#endif // OPENSWMM_ENGINE_FV_SECTION_GEOMETRY_HPP