LinkSubtypes.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_LINK_SUBTYPES_HPP
#define OPENSWMM_ENGINE_LINK_SUBTYPES_HPP
 
#include <algorithm>
#include <vector>
#include <cstdint>
#include <string>
#include <type_traits>
 
#include "LinkData.hpp"
 
namespace openswmm {
 
// ============================================================================
// ConduitData — side-table for CONDUIT links
// ============================================================================

struct ConduitData {
    std::vector<int> link_idx; 
 
    std::vector<double> roughness;
    std::vector<double> length; 
    std::vector<double> true_length;
    std::vector<double> slope;
    std::vector<double> mod_length;
    std::vector<int> barrels;
    std::vector<double> beta;
    std::vector<double> rough_factor;
    std::vector<double> q_full;
    std::vector<double> q_max;
    std::vector<double> loss_inlet;
    std::vector<double> loss_outlet;
    std::vector<double> loss_avg;
    std::vector<double> seep_rate;
    std::vector<int> culvert_code;
 
    // Phase 6 Stage B: mutable per-step conduit state moved off wide LinkData.
    std::vector<double> evap_loss_rate;       
    std::vector<double> seep_loss_rate;       
    std::vector<uint8_t> normal_flow_limited; 
    std::vector<uint8_t> inlet_control;       
    std::vector<int8_t> full_state;           
 
    // G-X4 (2026-09-20): the two-zone aquifer's frozen coupling for this
    // routing step, published by SurfaceRouter2D before routing and read by
    // the two loss laws (DynamicWave and routing::computeConduitLosses).
    // `gw_coupled == 0` on every conduit of every deck without a
    // `[2D_AQUIFER]` + `LINK_SEEPAGE TWO_WAY`, and the laws then take their
    // legacy branch unchanged.
    std::vector<uint8_t> gw_coupled; 
    std::vector<double>
        gw_head_rel; 
    std::vector<double> gw_kc; 
    std::vector<double>
        gw_dc; 
    std::vector<double>
        gw_gain_max; 

    int count() const noexcept { return static_cast<int>(link_idx.size()); }

    void clear() noexcept {
        link_idx.clear();
        roughness.clear();
        length.clear();
        true_length.clear();
        slope.clear();
        mod_length.clear();
        barrels.clear();
        beta.clear();
        rough_factor.clear();
        q_full.clear();
        q_max.clear();
        loss_inlet.clear();
        loss_outlet.clear();
        loss_avg.clear();
        seep_rate.clear();
        culvert_code.clear();
        evap_loss_rate.clear();
        seep_loss_rate.clear();
        normal_flow_limited.clear();
        inlet_control.clear();
        full_state.clear();
        gw_coupled.clear();
        gw_head_rel.clear();
        gw_kc.clear(); // G-X4
        gw_dc.clear();
        gw_gain_max.clear();
    }

    void reserve(int n) {
        const auto un = static_cast<std::size_t>(n);
        link_idx.reserve(un);
        roughness.reserve(un);
        length.reserve(un);
        true_length.reserve(un);
        slope.reserve(un);
        mod_length.reserve(un);
        barrels.reserve(un);
        beta.reserve(un);
        rough_factor.reserve(un);
        q_full.reserve(un);
        q_max.reserve(un);
        loss_inlet.reserve(un);
        loss_outlet.reserve(un);
        loss_avg.reserve(un);
        seep_rate.reserve(un);
        culvert_code.reserve(un);
        evap_loss_rate.reserve(un);
        seep_loss_rate.reserve(un);
        normal_flow_limited.reserve(un);
        inlet_control.reserve(un);
        full_state.reserve(un);
        gw_coupled.reserve(un);
        gw_head_rel.reserve(un);
        gw_kc.reserve(un); // G-X4
        gw_dc.reserve(un);
        gw_gain_max.reserve(un);
    }

    int add_default(int i) {
        const auto p = static_cast<std::ptrdiff_t>(
            std::lower_bound(link_idx.begin(), link_idx.end(), i) - link_idx.begin());
        link_idx.insert(link_idx.begin() + p, i);
        roughness.insert(roughness.begin() + p, 0.01);
        length.insert(length.begin() + p, 0.0);
        true_length.insert(true_length.begin() + p, 0.0);
        slope.insert(slope.begin() + p, 0.0);
        mod_length.insert(mod_length.begin() + p, 0.0);
        barrels.insert(barrels.begin() + p, 1);
        beta.insert(beta.begin() + p, 0.0);
        rough_factor.insert(rough_factor.begin() + p, 0.0);
        q_full.insert(q_full.begin() + p, 0.0);
        q_max.insert(q_max.begin() + p, 0.0);
        loss_inlet.insert(loss_inlet.begin() + p, 0.0);
        loss_outlet.insert(loss_outlet.begin() + p, 0.0);
        loss_avg.insert(loss_avg.begin() + p, 0.0);
        seep_rate.insert(seep_rate.begin() + p, 0.0);
        culvert_code.insert(culvert_code.begin() + p, 0);
        evap_loss_rate.insert(evap_loss_rate.begin() + p, 0.0);
        seep_loss_rate.insert(seep_loss_rate.begin() + p, 0.0);
        normal_flow_limited.insert(normal_flow_limited.begin() + p, uint8_t{0});
        inlet_control.insert(inlet_control.begin() + p, uint8_t{0});
        full_state.insert(full_state.begin() + p, int8_t{0});
        gw_coupled.insert(gw_coupled.begin() + p, uint8_t{0}); // G-X4
        gw_head_rel.insert(gw_head_rel.begin() + p, 0.0);
        gw_kc.insert(gw_kc.begin() + p, 0.0);
        gw_dc.insert(gw_dc.begin() + p, 0.0);
        gw_gain_max.insert(gw_gain_max.begin() + p, 0.0);
        return static_cast<int>(p);
    }

    void erase_at(int r) {
        const auto p = static_cast<std::ptrdiff_t>(r);
        link_idx.erase(link_idx.begin() + p);
        roughness.erase(roughness.begin() + p);
        length.erase(length.begin() + p);
        true_length.erase(true_length.begin() + p);
        slope.erase(slope.begin() + p);
        mod_length.erase(mod_length.begin() + p);
        barrels.erase(barrels.begin() + p);
        beta.erase(beta.begin() + p);
        rough_factor.erase(rough_factor.begin() + p);
        q_full.erase(q_full.begin() + p);
        q_max.erase(q_max.begin() + p);
        loss_inlet.erase(loss_inlet.begin() + p);
        loss_outlet.erase(loss_outlet.begin() + p);
        loss_avg.erase(loss_avg.begin() + p);
        seep_rate.erase(seep_rate.begin() + p);
        culvert_code.erase(culvert_code.begin() + p);
        evap_loss_rate.erase(evap_loss_rate.begin() + p);
        seep_loss_rate.erase(seep_loss_rate.begin() + p);
        normal_flow_limited.erase(normal_flow_limited.begin() + p);
        inlet_control.erase(inlet_control.begin() + p);
        full_state.erase(full_state.begin() + p);
    }
};
 
// ============================================================================
// PumpData — side-table for PUMP links
// ============================================================================
 
struct PumpData {
    std::vector<int> link_idx;
    std::vector<int> curve;          
    std::vector<uint8_t> init_state; 
    std::vector<double> startup;     
    std::vector<double> shutoff;     
    std::vector<int> curve_type;     

    int count() const noexcept { return static_cast<int>(link_idx.size()); }

    void clear() noexcept {
        link_idx.clear();
        curve.clear();
        init_state.clear();
        startup.clear();
        shutoff.clear();
        curve_type.clear();
    }

    void reserve(int n) {
        const auto un = static_cast<std::size_t>(n);
        link_idx.reserve(un);
        curve.reserve(un);
        init_state.reserve(un);
        startup.reserve(un);
        shutoff.reserve(un);
        curve_type.reserve(un);
    }

    int add_default(int i) {
        const auto p = static_cast<std::ptrdiff_t>(
            std::lower_bound(link_idx.begin(), link_idx.end(), i) - link_idx.begin());
        link_idx.insert(link_idx.begin() + p, i);
        curve.insert(curve.begin() + p, -1);
        init_state.insert(init_state.begin() + p, uint8_t{0});
        startup.insert(startup.begin() + p, 0.0);
        shutoff.insert(shutoff.begin() + p, 0.0);
        curve_type.insert(curve_type.begin() + p, -1);
        return static_cast<int>(p);
    }

    void erase_at(int r) {
        const auto p = static_cast<std::ptrdiff_t>(r);
        link_idx.erase(link_idx.begin() + p);
        curve.erase(curve.begin() + p);
        init_state.erase(init_state.begin() + p);
        startup.erase(startup.begin() + p);
        shutoff.erase(shutoff.begin() + p);
        curve_type.erase(curve_type.begin() + p);
    }
};
 
// ============================================================================
// OrificeData — side-table for ORIFICE links  (wide param1 → orifice_type)
// ============================================================================
 
struct OrificeData {
    std::vector<int> link_idx;
    std::vector<double> orifice_type; 
    std::vector<double> cd;           
    std::vector<double> orate;        
    int count() const noexcept { return static_cast<int>(link_idx.size()); }

    void clear() noexcept {
        link_idx.clear();
        orifice_type.clear();
        cd.clear();
        orate.clear();
    }

    void reserve(int n) {
        const auto un = static_cast<std::size_t>(n);
        link_idx.reserve(un);
        orifice_type.reserve(un);
        cd.reserve(un);
        orate.reserve(un);
    }

    int add_default(int i) {
        const auto p = static_cast<std::ptrdiff_t>(
            std::lower_bound(link_idx.begin(), link_idx.end(), i) - link_idx.begin());
        link_idx.insert(link_idx.begin() + p, i);
        orifice_type.insert(orifice_type.begin() + p, 0.0);
        cd.insert(cd.begin() + p, 0.0);
        orate.insert(orate.begin() + p, 0.0);
        return static_cast<int>(p);
    }

    void erase_at(int r) {
        const auto p = static_cast<std::ptrdiff_t>(r);
        link_idx.erase(link_idx.begin() + p);
        orifice_type.erase(orifice_type.begin() + p);
        cd.erase(cd.begin() + p);
        orate.erase(orate.begin() + p);
    }
};
 
// ============================================================================
// WeirData — side-table for WEIR links  (param1 → weir_type, param2 → end_con)
// ============================================================================
 
struct WeirData {
    std::vector<int> link_idx;
    std::vector<double>
        weir_type; 
    std::vector<double> cd;               
    std::vector<double> end_contractions; 
    std::vector<double> crest_height;
    std::vector<double> cd2; 
    std::vector<uint8_t>
        can_surcharge; 
    std::vector<double>
        road_width; 
    std::vector<int8_t>
        road_surface; 
    std::vector<int>
        cd_curve; 

    int count() const noexcept { return static_cast<int>(link_idx.size()); }

    void clear() noexcept {
        link_idx.clear();
        weir_type.clear();
        cd.clear();
        end_contractions.clear();
        crest_height.clear();
        cd2.clear();
        can_surcharge.clear();
        road_width.clear();
        road_surface.clear();
        cd_curve.clear();
    }

    void reserve(int n) {
        const auto un = static_cast<std::size_t>(n);
        link_idx.reserve(un);
        weir_type.reserve(un);
        cd.reserve(un);
        end_contractions.reserve(un);
        crest_height.reserve(un);
        cd2.reserve(un);
        can_surcharge.reserve(un);
        road_width.reserve(un);
        road_surface.reserve(un);
        cd_curve.reserve(un);
    }

    int add_default(int i) {
        const auto p = static_cast<std::ptrdiff_t>(
            std::lower_bound(link_idx.begin(), link_idx.end(), i) - link_idx.begin());
        link_idx.insert(link_idx.begin() + p, i);
        weir_type.insert(weir_type.begin() + p, 0.0);
        cd.insert(cd.begin() + p, 0.0);
        end_contractions.insert(end_contractions.begin() + p, 0.0);
        crest_height.insert(crest_height.begin() + p, 0.0);
        cd2.insert(cd2.begin() + p, 0.0);
        can_surcharge.insert(can_surcharge.begin() + p, uint8_t{1});
        road_width.insert(road_width.begin() + p, 0.0);
        road_surface.insert(road_surface.begin() + p, int8_t{0});
        cd_curve.insert(cd_curve.begin() + p, -1);
        return static_cast<int>(p);
    }

    void erase_at(int r) {
        const auto p = static_cast<std::ptrdiff_t>(r);
        link_idx.erase(link_idx.begin() + p);
        weir_type.erase(weir_type.begin() + p);
        cd.erase(cd.begin() + p);
        end_contractions.erase(end_contractions.begin() + p);
        crest_height.erase(crest_height.begin() + p);
        cd2.erase(cd2.begin() + p);
        can_surcharge.erase(can_surcharge.begin() + p);
        road_width.erase(road_width.begin() + p);
        road_surface.erase(road_surface.begin() + p);
        cd_curve.erase(cd_curve.begin() + p);
    }
};
 
// ============================================================================
// OutletData — side-table for OUTLET links
//   param1 → outlet_type (rating: 0 FUNC/HEAD,1 FUNC/DEPTH,2 TAB/HEAD,3 TAB/DEPTH)
//   cd → coeff (functional C1), param2 → expon (functional C2)
// ============================================================================
 
struct OutletData {
    std::vector<int> link_idx;
    std::vector<double> outlet_type;  
    std::vector<double> crest_height; 
    std::vector<double> coeff;        
    std::vector<double> expon;        
    std::vector<int>
        curve; 

    int count() const noexcept { return static_cast<int>(link_idx.size()); }

    void clear() noexcept {
        link_idx.clear();
        outlet_type.clear();
        crest_height.clear();
        coeff.clear();
        expon.clear();
        curve.clear();
    }

    void reserve(int n) {
        const auto un = static_cast<std::size_t>(n);
        link_idx.reserve(un);
        outlet_type.reserve(un);
        crest_height.reserve(un);
        coeff.reserve(un);
        expon.reserve(un);
        curve.reserve(un);
    }

    int add_default(int i) {
        const auto p = static_cast<std::ptrdiff_t>(
            std::lower_bound(link_idx.begin(), link_idx.end(), i) - link_idx.begin());
        link_idx.insert(link_idx.begin() + p, i);
        outlet_type.insert(outlet_type.begin() + p, 0.0);
        crest_height.insert(crest_height.begin() + p, 0.0);
        coeff.insert(coeff.begin() + p, 0.0);
        expon.insert(expon.begin() + p, 0.0);
        curve.insert(curve.begin() + p, -1);
        return static_cast<int>(p);
    }

    void erase_at(int r) {
        const auto p = static_cast<std::ptrdiff_t>(r);
        link_idx.erase(link_idx.begin() + p);
        outlet_type.erase(outlet_type.begin() + p);
        crest_height.erase(crest_height.begin() + p);
        coeff.erase(coeff.begin() + p);
        expon.erase(expon.begin() + p);
        curve.erase(curve.begin() + p);
    }
};
 
// ============================================================================
// LinkSubtypes — container + reverse map + incremental maintenance
// ============================================================================

struct LinkSubtypes {
    ConduitData conduits;
    PumpData pumps;
    OrificeData orifices;
    WeirData weirs;
    OutletData outlets;

    std::vector<int> subtype_row;

    void build_default_rows(const LinkData& L) {
        clear();
        const int n = L.count();
        for (int i = 0; i < n; ++i) {
            switch (L.type[static_cast<std::size_t>(i)]) {
            case LinkType::CONDUIT:
                conduits.add_default(i);
                break;
            case LinkType::PUMP:
                pumps.add_default(i);
                break;
            case LinkType::ORIFICE:
                orifices.add_default(i);
                break;
            case LinkType::WEIR:
                weirs.add_default(i);
                break;
            case LinkType::OUTLET:
                outlets.add_default(i);
                break;
            }
        }
        rebuild_index(n);
    }

    void clear() noexcept {
        conduits.clear();
        pumps.clear();
        orifices.clear();
        weirs.clear();
        outlets.clear();
        subtype_row.clear();
    }

    void rebuild_index(int n_links) {
        subtype_row.assign(static_cast<std::size_t>(n_links), -1);
        auto idx = [&](const std::vector<int>& keys) {
            for (int r = 0; r < static_cast<int>(keys.size()); ++r)
                subtype_row[static_cast<std::size_t>(keys[static_cast<std::size_t>(r)])] = r;
        };
        idx(conduits.link_idx);
        idx(pumps.link_idx);
        idx(orifices.link_idx);
        idx(weirs.link_idx);
        idx(outlets.link_idx);
    }

    int set_link_type(LinkData& links, int i, LinkType t) {
        const auto ui = static_cast<std::size_t>(i);
        if (i >= static_cast<int>(subtype_row.size()))
            subtype_row.resize(static_cast<std::size_t>(i) + 1, -1);
 
        const LinkType old = links.type[ui];
        const int existing = subtype_row[ui];
        if (old == t && existing >= 0)
            return existing; // already this subtype with a row — idempotent
 
        bool shifted = false;
        if (existing >= 0) { // remove the old subtype row (real re-type/convert)
            switch (old) {
            case LinkType::CONDUIT:
                conduits.erase_at(existing);
                break;
            case LinkType::PUMP:
                pumps.erase_at(existing);
                break;
            case LinkType::ORIFICE:
                orifices.erase_at(existing);
                break;
            case LinkType::WEIR:
                weirs.erase_at(existing);
                break;
            case LinkType::OUTLET:
                outlets.erase_at(existing);
                break;
            }
            shifted = true;
        }
 
        links.type[ui] = t;
 
        int row = -1;
        switch (t) {
        case LinkType::CONDUIT:
            row = conduits.add_default(i);
            if (row != conduits.count() - 1)
                shifted = true;
            break;
        case LinkType::PUMP:
            row = pumps.add_default(i);
            if (row != pumps.count() - 1)
                shifted = true;
            break;
        case LinkType::ORIFICE:
            row = orifices.add_default(i);
            if (row != orifices.count() - 1)
                shifted = true;
            break;
        case LinkType::WEIR:
            row = weirs.add_default(i);
            if (row != weirs.count() - 1)
                shifted = true;
            break;
        case LinkType::OUTLET:
            row = outlets.add_default(i);
            if (row != outlets.count() - 1)
                shifted = true;
            break;
        }
 
        if (shifted) {
            rebuild_index(links.count());
            return subtype_row[ui];
        }
        subtype_row[ui] = row; // O(1) ascending-parse / end-insert path
        return row;
    }

    void erase_link(int i, int n_after) {
        auto drop = [i](auto& tbl) -> bool {
            for (int r = 0; r < tbl.count(); ++r)
                if (tbl.link_idx[static_cast<std::size_t>(r)] == i) {
                    tbl.erase_at(r);
                    return true;
                }
            return false;
        };
        drop(conduits) || drop(pumps) || drop(orifices) || drop(weirs) || drop(outlets);
 
        auto shift = [i](std::vector<int>& keys) {
            for (auto& k : keys)
                if (k > i)
                    --k;
        };
        shift(conduits.link_idx);
        shift(pumps.link_idx);
        shift(orifices.link_idx);
        shift(weirs.link_idx);
        shift(outlets.link_idx);
 
        rebuild_index(n_after);
    }
 
    // O(1) row helpers (validate link_idx[r] == i, like NodeSubtypes::*_row).
    int conduit_row(int i) const noexcept { return row_in(i, conduits.link_idx); }
    int pump_row(int i) const noexcept { return row_in(i, pumps.link_idx); }
    int orifice_row(int i) const noexcept { return row_in(i, orifices.link_idx); }
    int weir_row(int i) const noexcept { return row_in(i, weirs.link_idx); }
    int outlet_row(int i) const noexcept { return row_in(i, outlets.link_idx); }
 
private:
    int row_in(int i, const std::vector<int>& keys) const noexcept {
        if (i >= 0 && i < static_cast<int>(subtype_row.size())) {
            const int r = subtype_row[static_cast<std::size_t>(i)];
            if (r >= 0 && r < static_cast<int>(keys.size()) &&
                keys[static_cast<std::size_t>(r)] == i)
                return r;
        }
        return -1;
    }
};
 
} // namespace openswmm
 
#endif // OPENSWMM_ENGINE_LINK_SUBTYPES_HPP