MeshData.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_2D_MESH_DATA_HPP
#define OPENSWMM_ENGINE_2D_MESH_DATA_HPP
 
#include <vector>
#include <string>
#include <cstdint>
 
namespace openswmm::twoD {

inline constexpr int kMaxCellVerts = 4;

struct MeshData {
 
    // -----------------------------------------------------------------------
    // Vertex arrays — indexed by vertex index [0, n_vertices)
    // -----------------------------------------------------------------------
 
    std::vector<double> vx;             
    std::vector<double> vy;             
    std::vector<double> vz;             
    std::vector<std::string> vtag;      
 
    // -----------------------------------------------------------------------
    // Cell static properties — indexed by cell index [0, n_triangles())
    // -----------------------------------------------------------------------

    std::vector<uint8_t> cell_nv;
    std::vector<int> cell_v;
    std::vector<int> cell_nbr;
 
    // Precomputed geometry (per cell)
    std::vector<double> tri_area;       
    std::vector<double> tri_cx;         
    std::vector<double> tri_cy;         
    std::vector<double> tri_cz;         
 
    // Edge geometry — flat 2D: [cell * kMaxCellVerts + k]
    std::vector<double> edge_length;    
    std::vector<double> edge_nx;        
    std::vector<double> edge_ny;        
    std::vector<double> edge_mx;        
    std::vector<double> edge_my;        
    std::vector<double> edge_mz;        
    std::vector<double> edge_dist_c;

    std::vector<double> quad_vfr_z;
    std::vector<double> quad_vfr_a;

    std::vector<double> edge_conveyance;
 
    // Surface properties
    std::vector<double> mannings_n;     
    std::vector<double> tri_init_depth; 
    std::vector<double> tri_init_u;     
    std::vector<double> tri_init_v;     
    std::vector<std::string> tri_tag;   
 
    // -----------------------------------------------------------------------
    // Vertex reconstruction stencil (pseudo-Laplacian weights)
    // -----------------------------------------------------------------------
    // Stored as CSR (compressed sparse row) for variable stencil sizes
    std::vector<int>    vert_stencil_ptr;   
    std::vector<int>    vert_stencil_idx;   
    std::vector<double> vert_stencil_wt;    
 
    // -----------------------------------------------------------------------
    // Coupling maps
    // -----------------------------------------------------------------------
 
    std::vector<int> vert_coupled_node;     
    std::vector<int> tri_coupled_node;      
 
    // Coupling parameters (per vertex/triangle coupling point)
    std::vector<double> vert_coupling_cd;       
    std::vector<double> vert_coupling_area;     
    std::vector<uint8_t> vert_coupling_area_set;
    std::vector<double> tri_coupling_cd;        
    std::vector<double> tri_coupling_area;      
 
    // Deferred resolution names (populated during parsing, cleared after resolve)
    std::vector<std::string> vert_coupled_node_name;
    std::vector<std::string> tri_coupled_node_name;

    struct TriCouplingRow {
        int         tri  = -1;    
        int         node = -1;    
        std::string node_name;    
        double      cd   = 0.65;  
        double      area = 1.0;   
        bool        area_set = false; 
    };
    std::vector<TriCouplingRow> tri_couplings;
 
    // -----------------------------------------------------------------------
    // Capacity queries
    // -----------------------------------------------------------------------
 
    int n_vertices()  const noexcept { return static_cast<int>(vx.size()); }
    int n_triangles() const noexcept { return static_cast<int>(cell_nv.size()); }
    int n_cells()     const noexcept { return static_cast<int>(cell_nv.size()); }
    int n_quads() const noexcept {
        int q = 0;
        for (uint8_t nv : cell_nv) q += (nv == 4);
        return q;
    }

    int edge_stride() const noexcept { return n_quads() > 0 ? kMaxCellVerts : 3; }
    int n_edge_slots() const noexcept { return n_cells() * kMaxCellVerts; }
 
    // -----------------------------------------------------------------------
    // Cell / edge-slot accessors
    // -----------------------------------------------------------------------

    static constexpr int slot(int c, int k) noexcept { return c * kMaxCellVerts + k; }
    static constexpr int slot_local(int s) noexcept { return s % kMaxCellVerts; }
    static constexpr int slot_cell(int s)  noexcept { return s / kMaxCellVerts; }
 
    int cell_vertex_count(int c) const noexcept { return cell_nv[static_cast<std::size_t>(c)]; }
    int cell_vertex(int c, int k) const noexcept { return cell_v[static_cast<std::size_t>(slot(c, k))]; }
    int cell_neighbour(int c, int k) const noexcept { return cell_nbr[static_cast<std::size_t>(slot(c, k))]; }
    bool is_quad(int c) const noexcept { return cell_nv[static_cast<std::size_t>(c)] == 4; }

    void cell_edge_vertices(int c, int k, int& a, int& b) const noexcept {
        const int nv = cell_nv[static_cast<std::size_t>(c)];
        a = cell_v[static_cast<std::size_t>(slot(c, (k + 1) % nv))];
        b = cell_v[static_cast<std::size_t>(slot(c, (k + 2) % nv))];
    }
 
    void set_triangle(int c, int v0, int v1, int v2) noexcept {
        cell_nv[static_cast<std::size_t>(c)] = 3;
        cell_v[static_cast<std::size_t>(slot(c, 0))] = v0;
        cell_v[static_cast<std::size_t>(slot(c, 1))] = v1;
        cell_v[static_cast<std::size_t>(slot(c, 2))] = v2;
        cell_v[static_cast<std::size_t>(slot(c, 3))] = -1;
    }
    void set_quad(int c, int v0, int v1, int v2, int v3) noexcept {
        cell_nv[static_cast<std::size_t>(c)] = 4;
        cell_v[static_cast<std::size_t>(slot(c, 0))] = v0;
        cell_v[static_cast<std::size_t>(slot(c, 1))] = v1;
        cell_v[static_cast<std::size_t>(slot(c, 2))] = v2;
        cell_v[static_cast<std::size_t>(slot(c, 3))] = v3;
    }
 
    // -----------------------------------------------------------------------
    // Resize / allocation
    // -----------------------------------------------------------------------
 
    void resize_vertices(int nv) {
        auto n = static_cast<std::size_t>(nv);
        vx.resize(n, 0.0);
        vy.resize(n, 0.0);
        vz.resize(n, 0.0);
        vtag.resize(n);
        vert_coupled_node.resize(n, -1);
        vert_coupled_node_name.resize(n);
        vert_coupling_cd.resize(n, 0.65);
        vert_coupling_area.resize(n, 1.0);
        vert_coupling_area_set.resize(n, 0);
    }

    void resize_triangles(int nt) {
        auto n = static_cast<std::size_t>(nt);
        const auto old_n = cell_nv.size();
        cell_nv.resize(n, 3);
        const auto n4 = n * static_cast<std::size_t>(kMaxCellVerts);
        cell_v.resize(n4, 0);
        cell_nbr.resize(n4, -1);
        // Padding slot of every NEW triangle: vertex −1, neighbour −2.
        for (std::size_t c = old_n; c < n; ++c) {
            cell_v[c * kMaxCellVerts + 3]   = -1;
            cell_nbr[c * kMaxCellVerts + 3] = -2;
        }
        tri_area.resize(n, 0.0);
        tri_cx.resize(n, 0.0);
        tri_cy.resize(n, 0.0);
        tri_cz.resize(n, 0.0);
 
        edge_length.resize(n4, 0.0);
        edge_nx.resize(n4, 0.0);
        edge_ny.resize(n4, 0.0);
        edge_mx.resize(n4, 0.0);
        edge_my.resize(n4, 0.0);
        edge_mz.resize(n4, 0.0);
        edge_dist_c.resize(n4, 0.0);
        edge_conveyance.resize(n4, 1.0);  // §11A — default unrestricted
 
        mannings_n.resize(n, 0.035);
        tri_init_depth.resize(n, 0.0);
        tri_init_u.resize(n, 0.0);
        tri_init_v.resize(n, 0.0);
        tri_tag.resize(n);
        tri_coupled_node.resize(n, -1);
        tri_coupled_node_name.resize(n);
        tri_coupling_cd.resize(n, 0.65);
        tri_coupling_area.resize(n, 1.0);
    }
};
 
} // namespace openswmm::twoD
 
#endif // OPENSWMM_ENGINE_2D_MESH_DATA_HPP