"""Три нативных эталона FEMM 5: электростатика, теплопроводность и ток.

Назначение: построить три линейные модели с известным одномерным решением,
решить каждую задачу и сравнить значения в центре с аналитическими числами.
Среда: встроенный Python FEMM 5 с доступными командами ei/eo, hi/ho, ci/co.
Внешние пакеты и pyFEMM не требуются. Геометрия задаётся в миллиметрах,
полевые результаты считываются в единицах СИ.

Вход: численные параметры непосредственно в разделах ниже. Выход: три модели,
соответствующие решения и result.json; поля real/imag сохраняют комплексные
значения без потери мнимой части. Проверяемые компоненты и их единицы указаны
рядом с expected. Допуск 10⁻⁶ относится только к этим линейным эталонам.

До запуска выберите отдельный пустой рабочий каталог FEMM: имена файлов
фиксированы и заданы относительно текущего каталога. Запуск создаёт новые
расчёты и может заменить одноимённые файлы. Сохраните открытые документы.
Успешное сравнение относится к введённым идеализированным моделям.
"""

import json, math, time

# Общий журнал накапливает результаты всех трёх задач и время выполнения, с.
out = {"cases": [], "passed": False}
started = time.time()

# 1. Электростатика: прямоугольник 10 × 10 мм, глубина 100 мм, ε_r = 2.
# newdocument(1) открывает электростатический препроцессор; planar задаёт
# плоское поперечное сечение. Напряжение слева 0 В, справа 100 В.
newdocument(1)
ei_probdef("millimeters", "planar", 1e-9, 100, 30)
ei_addmaterial("Dielectric", 2, 2, 0)
ei_addboundprop("Left", 0, 0, 0, 0, 0)
ei_addboundprop("Right", 100, 0, 0, 0, 0)

# Последовательно замыкаем контур. На верхней и нижней сторонах остаётся
# естественное условие нулевой нормальной электрической индукции.
ei_addnode(0, 0)
ei_addnode(10, 0)
ei_addnode(10, 10)
ei_addnode(0, 10)
ei_addsegment(0, 0, 10, 0)
ei_addsegment(10, 0, 10, 10)
ei_addsegment(10, 10, 0, 10)
ei_addsegment(0, 10, 0, 0)
ei_selectsegment(0, 5.0)
ei_setsegmentprop("Left", 0, 1, 0, 0, "<None>")
ei_clearselected()
ei_selectsegment(10, 5.0)
ei_setsegmentprop("Right", 0, 1, 0, 0, "<None>")
ei_clearselected()

# Метка назначает диэлектрик всей области и предельный размер элемента 2 мм.
ei_addblocklabel(5.0, 5.0)
ei_selectlabel(5.0, 5.0)
ei_setblockprop("Dielectric", 0, 2.0, 0)
ei_clearselected()

# Сначала сохраняем вход, затем решаем и загружаем результат для команд eo.
ei_saveas('./electrostatic.fee')
ei_analyze(1)
ei_loadsolution()
values = eo_getpointvalues(5.0, 5.0)

# Индексы Python начинаются с нуля: [0] — потенциал, В; [3] — E_x, В/м.
# В центре V = 50 В; E_x = -dV/dx = -100/0,01 = -10000 В/м.
expected = {0: 50, 3: -10000}
errors = {
    str(k): abs(complex(values[k]) - v) / max(abs(v), 1e-12)
    for k, v in expected.items()
}
out["cases"].append({
    "physics": 'electrostatic',
    "model": 'electrostatic.fee',
    "point": [5.0, 5.0],
    "values": [{"real": complex(v).real, "imag": complex(v).imag} for v in values],
    "expected": {str(k): v for k, v in expected.items()},
    "relative_errors": errors,
    "passed": max(errors.values()) < 1e-6,
})

# Закрываем постпроцессор, затем входной документ перед следующим типом задачи.
eo_close()
ei_close()

# 2. Теплопроводность: прямоугольник 10 × 10 мм, глубина 100 мм.
# k_x = k_y = 10 Вт/(м·К), объёмного источника нет. Последний аргумент
# материала — c_v = 1 МДж/(м³·К); в этой стационарной задаче он не влияет.
newdocument(2)
hi_probdef("millimeters", "planar", 1e-9, 100, 30)
hi_addmaterial("Solid", 10, 10, 0, 1)
hi_addboundprop("Left", 0, 373.15, 0, 0, 0, 0)
hi_addboundprop("Right", 0, 293.15, 0, 0, 0, 0)

# Температуры торцов заданы в кельвинах; верх и низ теплоизолированы.
hi_addnode(0, 0)
hi_addnode(10, 0)
hi_addnode(10, 10)
hi_addnode(0, 10)
hi_addsegment(0, 0, 10, 0)
hi_addsegment(10, 0, 10, 10)
hi_addsegment(10, 10, 0, 10)
hi_addsegment(0, 10, 0, 0)
hi_selectsegment(0, 5.0)
hi_setsegmentprop("Left", 0, 1, 0, 0, "<None>")
hi_clearselected()
hi_selectsegment(10, 5.0)
hi_setsegmentprop("Right", 0, 1, 0, 0, "<None>")
hi_clearselected()
hi_addblocklabel(5.0, 5.0)
hi_selectlabel(5.0, 5.0)
hi_setblockprop("Solid", 0, 2.0, 0)
hi_clearselected()
hi_saveas('./heat.feh')
hi_analyze(1)
hi_loadsolution()
values = ho_getpointvalues(5.0, 5.0)

# [0] — температура, К; [1] — плотность теплового потока F_x, Вт/м².
# T_mid = (373,15 + 293,15)/2; F_x = -k dT/dx = 80000 Вт/м².
expected = {0: 333.15, 1: 80000}
errors = {
    str(k): abs(complex(values[k]) - v) / max(abs(v), 1e-12)
    for k, v in expected.items()
}
out["cases"].append({
    "physics": 'heat',
    "model": 'heat.feh',
    "point": [5.0, 5.0],
    "values": [{"real": complex(v).real, "imag": complex(v).imag} for v in values],
    "expected": {str(k): v for k, v in expected.items()},
    "relative_errors": errors,
    "passed": max(errors.values()) < 1e-6,
})
ho_close()
hi_close()

# 3. Постоянный ток: проводник 100 × 1 мм, глубина 10 мм; частота 0 Гц.
# σ_x = σ_y = 58·10⁶ См/м, потенциал слева 0,1 В, справа 0 В.
newdocument(3)
ci_probdef("millimeters", "planar", 0, 1e-9, 10, 30)
ci_addmaterial("Copper", 58000000, 58000000, 1, 1, 0, 0)
ci_addboundprop("Left", 0.1, 0, 0, 0, 0)
ci_addboundprop("Right", 0, 0, 0, 0, 0)

# Замкнутый контур и отсутствие боковых источников дают одномерное поле.
ci_addnode(0, 0)
ci_addnode(100, 0)
ci_addnode(100, 1)
ci_addnode(0, 1)
ci_addsegment(0, 0, 100, 0)
ci_addsegment(100, 0, 100, 1)
ci_addsegment(100, 1, 0, 1)
ci_addsegment(0, 1, 0, 0)
ci_selectsegment(0, 0.5)
ci_setsegmentprop("Left", 0, 1, 0, 0, "<None>")
ci_clearselected()
ci_selectsegment(100, 0.5)
ci_setsegmentprop("Right", 0, 1, 0, 0, "<None>")
ci_clearselected()
ci_addblocklabel(50.0, 0.5)
ci_selectlabel(50.0, 0.5)
ci_setblockprop("Copper", 0, 0.2, 0)
ci_clearselected()
ci_saveas('./current.fec')
ci_analyze(1)
ci_loadsolution()
values = co_getpointvalues(50.0, 0.5)

# [0] — потенциал, В; [1] — плотность тока J_x, А/м², а не полный ток.
# E_x = 0,1/0,1 = 1 В/м; J_x = σE_x = 58·10⁶ А/м².
expected = {0: 0.05, 1: 58000000}
errors = {
    str(k): abs(complex(values[k]) - v) / max(abs(v), 1e-12)
    for k, v in expected.items()
}
out["cases"].append({
    "physics": 'current',
    "model": 'current.fec',
    "point": [50.0, 0.5],
    "values": [{"real": complex(v).real, "imag": complex(v).imag} for v in values],
    "expected": {str(k): v for k, v in expected.items()},
    "relative_errors": errors,
    "passed": max(errors.values()) < 1e-6,
})
co_close()
ci_close()

# Итог положителен, только если прошли все три сравнения. Ненулевые эталоны
# нормируют ошибки; max(..., 10⁻¹²) защищает деление при расширении набора.
out["passed"] = all(v["passed"] for v in out["cases"])
out["seconds"] = time.time() - started
with open('./result.json', "w") as f:
    json.dump(out, f)
print(json.dumps(out))
assert out["passed"]