Building/assets/code/hu_from_slide.py

"""Средний коэффициент Hu порового давления по моделям Slide2.

В ПК Prorock гидростатическое воздействие задаётся формулой P = gamma_w h H_u
(раздел «Учёт порового давления»), где h -- вертикальное расстояние от узла
до заданной поверхности грунтовых вод, H_u -- коэффициент вклада поровой
воды. Поровое давление, рассчитанное в Slide2 фильтрационным анализом
(метод Groundwater), заменяется одним коэффициентом H_u, наилучшим образом
воспроизводящим рассчитанное поле. Поверхность грунтовых вод в Prorock при
этом задаётся по рельефу дневной поверхности, поэтому h -- глубина узла
ниже рельефа.

Файл .slim -- zip-архив модели Slide2, содержащий три текстовых файла:

  .sli -- модель: удельный вес воды gamma_w (ключ gammaw), координаты вершин
          геометрии (секция vertices) и внешний контур (exterior), верхняя
          огибающая которого задаёт рельеф дневной поверхности;
  .slw -- сетка фильтрационного расчёта: координаты узлов (секция nodes);
  .w01 -- результаты расчёта: строки «номер узла, полный напор, пьезометри-
          ческий напор» в метрах; поровое давление u = gamma_w h_p, кПа.

Для каждого узла сетки вычисляются глубина h ниже рельефа и поровое
давление u, после чего строятся три оценки коэффициента:

  1. H_u по МНК через начало координат u = H_u gamma_w h (узлы с h >= 1 м):
     H_u = sum(u h) / (gamma_w sum(h^2)); минимизирует среднеквадратичную
     невязку поля давлений -- основное значение для переноса в Prorock;
  2. среднее и медиана поточечных отношений H_u,i = u / (gamma_w h) по узлам
     с h >= 5 м (порог отсекает неустойчивые отношения у поверхности);
  3. среднее H_u,i в полосе изобары 1 МПа (900--1100 кПа) -- аналог прежнего
     ручного способа, в котором коэффициент определялся по глубине залегания
     изобары 1000 кПа.

Запуск (из корня репозитория):
  python assets/code/hu_from_slide.py "E:\\...\\Исходники_Slide" [ещё пути...]
Аргументы командной строки -- файлы .slim или каталоги; для каталога
обрабатываются все файлы .slim в его подкаталогах. Результаты печатаются
в консоль и записываются в assets/tables/hu_from_slide.csv. Имя модели
«35 1.slim» разбирается как контур 2035 года и профильная линия 1.
Используется только стандартная библиотека.
"""
import csv
import math
import re
import statistics
import sys
import zipfile
from pathlib import Path

# Пороги глубины, м: для МНК и для поточечных отношений H_u,i.
H_MIN_FIT, H_MIN_POINT = 1.0, 5.0
# Полоса изобары 1 МПа для сопоставления с прежним ручным способом, кПа.
ISOBAR_LO, ISOBAR_HI = 900.0, 1100.0
# Префикс номера контура в имени файла -- год контура.
CONTOUR_YEARS = {"26": 2026, "31": 2031, "33": 2033, "35": 2035}

_NUM = r"[-+]?\d*\.?\d+(?:[eE][-+]?\d+)?"


def section(text, name):
    """Строки секции name. Заголовок -- слово с двоеточием в начале строки,
    данные -- всё до следующего заголовка (в .sli с отступом, в .slw нет)."""
    lines, inside = [], False
    for line in text.splitlines():
        if line[:1].isspace() or line[:1].isdigit():
            if inside:
                lines.append(line)
        else:
            inside = line.split(":")[0].strip() == name
    return lines


def read_vertices(sli):
    """Координаты вершин геометрии: {номер: (x, y)}."""
    rx = re.compile(rf"(\d+)\s+x:\s*({_NUM})\s+y:\s*({_NUM})")
    return {int(m[1]): (float(m[2]), float(m[3]))
            for m in map(rx.search, section(sli, "vertices")) if m}


def read_exterior(sli):
    """Номера вершин внешнего контура (замкнутого многоугольника)."""
    for line in section(sli, "exterior"):
        m = re.search(r"\[([0-9,\s]+)\]", line)
        if m:
            return [int(v) for v in m[1].split(",")]
    return []


def surface_function(ids, verts):
    """Рельеф дневной поверхности: верхняя огибающая внешнего контура
    как функция от x."""
    ring = [verts[i] for i in ids if i in verts]
    ring += ring[:1]                      # замыкание контура
    segments = [(a, b) for a, b in zip(ring, ring[1:]) if a[0] != b[0]]

    def y_at(x):
        top = None
        for (x1, y1), (x2, y2) in segments:
            lo, hi = (x1, x2) if x1 < x2 else (x2, x1)
            if lo <= x <= hi:
                y = y1 + (y2 - y1)*(x - x1)/(x2 - x1)
                top = y if top is None or y > top else top
        return top
    return y_at


def read_nodes(slw):
    """Координаты узлов сетки фильтрационного расчёта: {номер: (x, y)}."""
    rx = re.compile(rf"(\d+)\s+x:\s*({_NUM})\s+y:\s*({_NUM})")
    return {int(m[1]): (float(m[2]), float(m[3]))
            for m in map(rx.search, section(slw, "nodes")) if m}


def read_heads(w01):
    """Напоры из результатов: {номер узла: (полный, пьезометрический)}, м."""
    heads = {}
    rx = re.compile(rf"(\d+)\s+({_NUM})\s+({_NUM})$")
    for line in w01.split("# Fluxes")[0].splitlines():  # блок напоров
        m = rx.match(line.strip())
        if m:
            heads[int(m[1])] = (float(m[2]), float(m[3]))
    return heads


def model_stats(path):
    """Оценки H_u по одной модели .slim."""
    with zipfile.ZipFile(path) as z:
        def pick(suffix):
            names = [n for n in z.namelist() if n.lower().endswith(suffix)]
            return z.read(names[0]).decode("cp1251", "replace") if names else ""
        sli, slw, w01 = pick(".sli"), pick(".slw"), pick(".w01")
    if not (sli and slw and w01):
        raise ValueError("нет .sli/.slw/.w01 -- модель без результатов расчёта")

    m = re.search(rf"gammaw:\s*({_NUM})", sli)
    gamma_w = float(m[1]) if m else 9.81

    verts = read_vertices(sli)
    y_surface = surface_function(read_exterior(sli), verts)
    nodes, heads = read_nodes(slw), read_heads(w01)

    # Глубина ниже рельефа и поровое давление каждого узла, м и кПа.
    # Контроль сопоставления узлов: полный напор минус высота узла равен
    # пьезометрическому напору (допускается невязка до 0,5 м).
    pts = []
    for n, (x, y) in nodes.items():
        if n in heads:
            total, ph = heads[n]
            if abs(total - y - ph) > 0.5:
                raise ValueError("напоры не сходятся с координатами узлов")
            y_top = y_surface(x)
            if y_top is not None and y_top - y >= 0:
                pts.append((y_top - y, gamma_w*ph))

    fit = [(h, u) for h, u in pts if h >= H_MIN_FIT]
    hu_lsm = sum(u*h for h, u in fit)/(gamma_w*sum(h*h for h, _ in fit))
    rms = math.sqrt(sum((u - hu_lsm*gamma_w*h)**2 for h, u in fit)/len(fit))

    point = [u/(gamma_w*h) for h, u in pts if h >= H_MIN_POINT and u > 0]
    isobar = [u/(gamma_w*h) for h, u in pts
              if h >= H_MIN_POINT and ISOBAR_LO <= u <= ISOBAR_HI]
    return {
        "n": len(fit), "h_max": max(h for h, _ in fit),
        "u_max": max(u for _, u in fit), "hu_lsm": hu_lsm, "rms": rms,
        "hu_mean": statistics.mean(point) if point else None,
        "hu_median": statistics.median(point) if point else None,
        "hu_isobar": statistics.mean(isobar) if isobar else None,
    }


def model_name(stem):
    """«35 1» -> (год контура, номер линии) по имени файла модели."""
    parts = stem.split()
    year = str(CONTOUR_YEARS.get(parts[0], parts[0]))
    return year, parts[1] if len(parts) > 1 else "--"


def fmt(value, digits=3):
    """Число с десятичной запятой и типографским знаком минус."""
    if value is None:
        return "--"
    return f"{value:.{digits}f}".replace(".", ",").replace("-", "\u2212")


def main():
    if len(sys.argv) < 2:
        sys.exit("укажите файлы .slim или каталоги с моделями Slide2")
    paths = []
    for arg in sys.argv[1:]:
        p = Path(arg)
        paths.extend(sorted(p.rglob("*.slim")) if p.is_dir() else [p])

    root = Path(__file__).resolve().parents[2]
    header = ("Год контура", "№ профильной линии", "Hu (МНК)",
              "Hu среднее", "Hu медиана", "Hu у изобары 1 МПа",
              "Узлов", "h_max, м", "u_max, кПа", "СКО поля, кПа")
    rows, values = [header], []
    print(f"{'год':>4} {'лин':>4} {'Hu МНК':>8} {'среднее':>8} "
          f"{'медиана':>8} {'1 МПа':>8} {'узлов':>6} "
          f"{'h_max':>6} {'u_max':>7} {'СКО':>7}")
    for path in paths:
        try:
            st = model_stats(path)
        except (ValueError, KeyError, zipfile.BadZipFile) as e:
            print(f"{path.stem}: пропущена ({e})")
            continue
        year, line = model_name(path.stem)
        print(f"{year:>4} {line:>4} {fmt(st['hu_lsm']):>8} "
              f"{fmt(st['hu_mean']):>8} {fmt(st['hu_median']):>8} "
              f"{fmt(st['hu_isobar']):>8} {st['n']:>6} "
              f"{fmt(st['h_max'], 0):>6} {fmt(st['u_max'], 0):>7} "
              f"{fmt(st['rms'], 0):>7}")
        rows.append((year, line, fmt(st["hu_lsm"]), fmt(st["hu_mean"]),
                     fmt(st["hu_median"]), fmt(st["hu_isobar"]), st["n"],
                     fmt(st["h_max"], 1), fmt(st["u_max"], 0),
                     fmt(st["rms"], 0)))
        values.append(st["hu_lsm"])

    if values:
        print(f"среднее Hu (МНК) по всем моделям: {fmt(statistics.mean(values))}")
        out = root/"assets/tables/hu_from_slide.csv"
        with open(out, "w", encoding="utf-8", newline="") as f:
            csv.writer(f, delimiter=";").writerows(rows)
        print(f"таблица записана: {out}")


if __name__ == "__main__":
    main()