Скругление электрода: интегральная и локальная сходимость

"""Скругление электрода и локальная напряжённость.
Среда: встроенный MicroPython FEMM 5; команды FEMM доступны глобально.
Геометрические размеры: миллиметры; поля и интегралы возвращаются в СИ.
Все материалы и проводники имеют русские имена. Каждая модель сохраняется,
закрывается и повторно открывается до нового расчёта, включая кириллицу.
Запуск создаёт новую папку calculated-NN рядом с программой; опубликованные
модели и предыдущие расчёты не перезаписываются. Все необходимые
исходные данные включены. Условия задачи и ограничения описаны в README.md.
"""
import os
import sys
import math
import json

saved = globals().get('__file__', '') or (sys.argv[0] if sys.argv else '')
base = os.path.dirname(os.path.abspath(saved)) if saved else os.getcwd()
index = 1
while os.path.exists(base + '/calculated-%02d' % index):
    index += 1
directory = base + '/calculated-%02d' % index
os.mkdir(directory)
records = []
checks = []
EPS0 = 8.85418781762e-12


def real(value):
    if isinstance(value, (list, tuple)):
        return real(value[0])
    return complex(value).real


def check(name, value, limit):
    checks.append({'name': name, 'value': value, 'limit': limit,
                   'passed': math.isfinite(value) and value <= limit})


def label(x, y, material, mesh, group=1):
    ei_addblocklabel(x, y)
    ei_selectlabel(x, y)
    ei_setblockprop(material, 0, mesh, group)
    ei_clearselected()


def circle(x, radius, conductor, angle):
    ei_addnode(x - radius, 0)
    ei_addnode(x + radius, 0)
    ei_addarc(x - radius, 0, x + radius, 0, 180, angle)
    ei_addarc(x + radius, 0, x - radius, 0, 180, angle)
    ei_selectarcsegment(x, radius)
    ei_selectarcsegment(x, -radius)
    ei_setarcsegmentprop(angle, '<None>', 0, 0, conductor)
    ei_clearselected()


def solve(model):
    ei_saveas(model)
    ei_close()
    opendocument(model)
    ei_analyze(1)
    ei_loadsolution()
    eo_smooth('off')


def electric_result(mesh, model, input_name, output_name, exact_C):
    voltage_in, charge_in = eo_getconductorproperties(input_name)
    voltage_out, charge_out = eo_getconductorproperties(output_name)
    voltage = real(voltage_in) - real(voltage_out)
    charge = real(charge_in)
    eo_groupselectblock(1)
    energy = real(eo_blockintegral(0))
    eo_clearblock()
    capacitance = charge / voltage
    energy_C = 2 * energy / voltage ** 2
    row = {'mesh_mm': mesh, 'model': os.path.basename(model),
           'nodes': eo_numnodes(), 'elements': eo_numelements(),
           'voltage_V': voltage, 'charge_C': charge, 'energy_J': energy,
           'capacitance_F': capacitance, 'energy_capacitance_F': energy_C,
           'analytic_capacitance_F': exact_C,
           'analytic_C_relative_error': abs(capacitance / exact_C - 1),
           'energy_charge_relative_error': abs(energy_C / capacitance - 1)}
    records.append(row)
    check('Энергия и заряд: ' + os.path.basename(model),
          row['energy_charge_relative_error'], 1e-5)
    return row


def circle_xy(cx, cy, radius, conductor, angle):
    ei_addnode(cx - radius, cy)
    ei_addnode(cx + radius, cy)
    ei_addarc(cx - radius, cy, cx + radius, cy, 180, angle)
    ei_addarc(cx + radius, cy, cx - radius, cy, 180, angle)
    ei_selectarcsegment(cx, cy + radius)
    ei_selectarcsegment(cx, cy - radius)
    ei_setarcsegmentprop(angle, '<None>', 0, 0, conductor)
    ei_clearselected()


def polygon(points, conductor='<None>', mesh=0):
    for x, y in points:
        ei_addnode(x, y)
    for k, (x, y) in enumerate(points):
        x2, y2 = points[(k + 1) % len(points)]
        ei_addsegment(x, y, x2, y2)
        if conductor != '<None>':
            ei_selectsegment((x + x2) / 2, (y + y2) / 2)
            ei_setsegmentprop('<None>', mesh, 0, 0, 0, conductor)
            ei_clearselected()


def terminal(x, y, name, mesh):
    ei_selectsegment(x, y)
    ei_setsegmentprop('<None>', mesh, 0, 0, 0, name)
    ei_clearselected()


def summary(mesh, model, source, ground, exact_C=None):
    v1, q1 = eo_getconductorproperties(source)
    v0, q0 = eo_getconductorproperties(ground)
    v1, v0, q1, q0 = real(v1), real(v0), real(q1), real(q0)
    voltage = v1 - v0
    eo_groupselectblock(1)
    energy = real(eo_blockintegral(0))
    eo_clearblock()
    capacitance = q1 / voltage
    row = {'mesh_mm': mesh, 'model': os.path.basename(model),
           'nodes': eo_numnodes(), 'elements': eo_numelements(),
           'voltage_V': voltage, 'charge_C': q1, 'energy_J': energy,
           'capacitance_F': capacitance,
           'energy_capacitance_F': 2 * energy / voltage ** 2}
    if exact_C is not None:
        row['analytic_capacitance_F'] = exact_C
        row['analytic_C_relative_error'] = abs(capacitance / exact_C - 1)
    records.append(row)
    check('Энергия и терминальный заряд: ' + os.path.basename(model),
          abs(2 * energy / voltage ** 2 / capacitance - 1), 1e-5)
    return row


def field(x, y):
    values = eo_getpointvalues(x, y)
    return math.sqrt(real(values[3]) ** 2 + real(values[4]) ** 2)


def rectangle_field(mesh, filename, epsx=2, epsy=2, voltage=100, width=10, height=10, depth=100):
    newdocument(1)
    ei_probdef('millimeters', 'planar', 1e-10, depth, 30)
    ei_addmaterial('Диэлектрик учебного образца', epsx, epsy, 0)
    ei_addconductorprop('Питающий электрод', voltage, 0, 1)
    ei_addconductorprop('Заземлённый электрод', 0, 0, 1)
    polygon(((0, 0), (width, 0), (width, height), (0, height)))
    terminal(0, height / 2, 'Питающий электрод', mesh)
    terminal(width, height / 2, 'Заземлённый электрод', mesh)
    label(width / 2, height / 2, 'Диэлектрик учебного образца', mesh)
    model = directory + '/' + filename + '.fee'
    solve(model)
    return model


def rounded_square(radius, angle, mesh):
    b = 3
    points = ((b - radius, b), (-b + radius, b), (-b, b - radius),
              (-b, -b + radius), (-b + radius, -b), (b - radius, -b),
              (b, -b + radius), (b, b - radius))
    for x, y in points:
        ei_addnode(x, y)
    for k in range(8):
        x, y = points[k]
        x2, y2 = points[(k + 1) % 8]
        if k % 2 == 0:
            ei_addsegment(x, y, x2, y2)
            terminal((x + x2) / 2, (y + y2) / 2, 'Скруглённый электрод', mesh)
        else:
            ei_addarc(x, y, x2, y2, 90, angle)
    p = b - radius + radius / math.sqrt(2)
    for x, y in ((p, p), (-p, p), (-p, -p), (p, -p)):
        ei_selectarcsegment(x, y)
    ei_setarcsegmentprop(angle, '<None>', 0, 0, 'Скруглённый электрод')
    ei_clearselected()


for radius in (0.2, 0.5, 1.0):
    for level, mesh in enumerate((0.15, 0.075, 0.0375)):
        newdocument(1)
        ei_probdef('millimeters', 'planar', 1e-10, 10, 30)
        ei_addmaterial('Воздух вокруг электрода', 1, 1, 0)
        ei_addconductorprop('Скруглённый электрод', 100, 0, 1)
        ei_addconductorprop('Заземлённый круговой экран', 0, 0, 1)
        angle = 4 / (2 ** level)
        rounded_square(radius, angle, mesh)
        circle(0, 12, 'Заземлённый круговой экран', angle)
        label(0, 0, '<No Mesh>', mesh, 0)
        label(6, 0, 'Воздух вокруг электрода', mesh)
        model = directory + '/rounding-%02d-%02d.fee' % (int(radius * 10), level + 1)
        solve(model)
        row = summary(mesh, model, 'Скруглённый электрод', 'Заземлённый круговой экран')
        sample = 3 - radius + (radius + 0.2) / math.sqrt(2)
        row['rounding_radius_mm'] = radius
        # Для локальной производной используем восстановленное узловое поле;
        # заряд и энергия выше получены независимо от этого сглаживания.
        eo_smooth('on')
        row['E_corner_sample_V_m'] = field(sample, sample)
        eo_close()
        ei_close()
    check('Сходимость ёмкости при радиусе ' + str(radius), abs(records[-1]['capacitance_F'] / records[-2]['capacitance_F'] - 1), 0.005)
    check('Сходимость поля в 0,2 мм от скругления ' + str(radius), abs(records[-1]['E_corner_sample_V_m'] / records[-2]['E_corner_sample_V_m'] - 1), 0.05)
check('Скругление снижает поле в заданной контрольной точке',
      max(0, records[8]['E_corner_sample_V_m'] / records[2]['E_corner_sample_V_m'] - 1), 0)

columns = {
    'mesh_mm': 'Размер сетки, мм', 'model': 'Файл модели',
    'nodes': 'Узлов', 'elements': 'Элементов',
    'voltage_V': 'Разность потенциалов, В', 'charge_C': 'Заряд, Кл',
    'energy_J': 'Энергия, Дж', 'capacitance_F': 'Ёмкость по заряду, Ф',
    'energy_capacitance_F': 'Ёмкость по энергии, Ф',
    'analytic_capacitance_F': 'Аналитическая ёмкость, Ф',
    'analytic_C_relative_error': 'Относительная ошибка ёмкости',
    'energy_charge_relative_error': 'Невязка энергии и заряда',
    'E_7p5mm_V_m': 'Поле при r=7,5 мм, В/м',
    'analytic_E_7p5mm_V_m': 'Аналитическое поле при r=7,5 мм, В/м',
    'E_150mm_V_m': 'Поле при r=150 мм, В/м',
    'analytic_E_150mm_V_m': 'Аналитическое поле при r=150 мм, В/м',
    'outer_radius_mm': 'Радиус внешней границы, мм',
    'insulation': 'Изоляция', 'analytic_Emax_kV_mm': 'Аналитический максимум поля, кВ/мм',
    'Emax_from_charge_kV_mm': 'Максимум поля по заряду, кВ/мм',
    'E_10mm_kV_mm': 'Поле при r=10 мм, кВ/мм',
    'E_25mm_kV_mm': 'Поле при r=25 мм, кВ/мм',
    'conductivity_S_m': 'Эквивалентная проводимость, См/м',
    'I_real_A': 'Активная компонента амплитуды тока, А',
    'I_imag_A': 'Квадратурная компонента амплитуды тока, А',
    'analytic_I_real_A': 'Аналитическая активная компонента, А',
    'analytic_I_imag_A': 'Аналитическая квадратурная компонента, А',
    'losses_W': 'Диэлектрические потери, Вт',
    'analytic_losses_W': 'Аналитические потери, Вт',
    'terminal_power_W': 'Средняя мощность электрода, Вт',
    'apparent_power_VA': 'Полная мощность, ВА',
    'U_RMS_V': 'Действующее напряжение, В', 'U_amplitude_V': 'Амплитуда напряжения, В'
}
columns.update({'orientation': 'Ориентация главных осей', 'E_center_V_m': 'Поле в центре, В/м', 'E_first_V_m': 'Поле первого слоя, В/м', 'E_second_V_m': 'Поле второго слоя, В/м', 'interface_V': 'Потенциал интерфейса, В', 'floating_voltage_V': 'Потенциал изолированного электрода, В', 'floating_charge_C': 'Заряд изолированного электрода, Кл', 'overlap_mm': 'Перекрытие электродов, мм', 'rounding_radius_mm': 'Радиус скругления, мм', 'E_corner_sample_V_m': 'Поле в 0,2 мм от скругления, В/м', 'screen': 'Исполнение', 'E_inner_sample_V_m': 'Поле в контрольной точке, В/м', 'excited_phase': 'Возбуждённая фаза', 'C_Aj_F': 'Коэффициент строки A, Ф', 'C_Bj_F': 'Коэффициент строки B, Ф', 'C_Cj_F': 'Коэффициент строки C, Ф'})
result = {'runtime': 'FEMM 5.0; встроенный MicroPython', 'records': records,
          'checks': checks, 'columns': columns,
          'passed': all(item['passed'] for item in checks)}
with open(directory + '/result.json', 'w') as output:
    json.dump(result, output)
print(json.dumps(result))
assert result['passed'], 'Не выполнен численный критерий'