Analiza parametryczna za pomocą Dlubal API

"""ParaS-VMStress_Thick-min.py

Powiązany model RFEM 6: "Para_Plate_Joint.rf6"
Model: https://www.dlubal.com/en/downloads-and-information/examples-and-tutorials/models-to-download/006159

Przeprowadzenie analizy parametrycznej za pomocą API Dlubal RFEM, zebranie metryk wyników, eksport i wykreślenie wyników.
oraz wykres.

Sekcje:
- Importy
- Konfiguracja
- Funkcje
- Wykonanie główne
- Wykresy
"""
# --------------------- Importy ---------------------
import itertools
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from dlubal.api import rfem


# --------------------- Konfiguracja ---------------------
parameter_ranges = {
    't_1': {'min': 8, 'max': 12, 'step': 2},
    't_2': {'min': 4, 'max': 8, 'step': 2},
} # grubości w mm
f_y = 235  # granica plastyczności w MPa

# --------------------- Funkcje ---------------------
def build_parameter_grid(param_ranges):
    """Buduje kompletny zbiór mutacji z wartości min/max/krok w mm."""
    columns = []
    value_lists = []
    for name, bounds in param_ranges.items():
        min_val = bounds['min']
        max_val = bounds['max']
        step = bounds['step']
        if step <= 0:
            raise ValueError(f"Szerokość kroku dla '{name}' musi być dodatnia.")
        if max_val < min_val:
            raise ValueError(f"Wartość maksymalna dla '{name}' musi być >= wartości minimalnej.")
        values = np.arange(min_val, max_val + step * 0.5, step)
        columns.append(name)
        value_lists.append(values)
    combinations = list(itertools.product(*value_lists))
    return pd.DataFrame(combinations, columns=columns)
       
def set_glpa(dl_app, p_name, p_value,):
    """Ustawia parametr globalny p_name na określoną wartość p_value"""
    params = dl_app.get_object_list([rfem.global_parameters.GlobalParameter()])
    for p in params:
        if p.name == p_name:
            p.value = p_value
            dl_app.update_object(p)
            check=True
            return True
    if not check:
        raise ValueError(f"  Błąd: Nie znaleziono parametru '{p_name}'!")

# --------------------- Wykonanie główne ---------------------
# Budowanie zbioru parametrów
para_dfo = build_parameter_grid(parameter_ranges)
para_dfo = para_dfo[para_dfo['t_1'] >= para_dfo['t_2']]
para_df = para_dfo / 1000  # mm -> m
print(f"Zbiór parametrów:\n {para_dfo}")

# Połącz z RFEM z bieżącym modelem i przeprowadź analizę parametryczną
with rfem.Application() as rf_app:
    # Sprawdź połączenie i wyświetl informacje o modelu
    app_info = rf_app.get_application_info()
    print("Informacje o aplikacji:", app_info)
    res_paras=pd.DataFrame(columns=['sigvm','sfz'])
    for i in para_df.index:
            # ustaw parametry globalne
            print(f"Ustaw mutację {i} z parametrami {[x for x in para_df.loc[i]]}")
            for p in para_df.columns:
                # print(f"Ustaw mutację {i} z parametrem {p} na {para_df.loc[i, p]}")
                if not set_glpa(rf_app, p, para_df.loc[i, p]):
                    break
            # Ponowna generacja siatki
            rf_app.delete_mesh()
            rf_app.generate_mesh(skip_warnings=True)
            # Uruchom obliczenia
            calculation = rf_app.calculate_all(skip_warnings=True)
            # Wyniki
            if calculation.succeeded:
                # Naprężenia von Misesa
                res_sigvm = rf_app.get_results(
                    results_type=rfem.results.STATIC_ANALYSIS_SURFACES_EQUIVALENT_STRESSES_MISES_MESH_NODES,
                    filters=[
                        rfem.results.ResultsFilter(column_id='loading', filter_expression='DS1'),
                        ]      
                    ).data
                sigvm = res_sigvm['sigma_eqv_mises'].max()/10**6 # N/m^2 -> MPa                
                # ciężar własny
                sfz = rf_app.get_result_table(
                    table=rfem.results.ResultTable.STATIC_ANALYSIS_SUMMARY_TABLE,
                    loading=rfem.ObjectId(no=1,object_type=rfem.OBJECT_TYPE_LOAD_CASE)
                ).data
                sfz = float(sfz.loc[6].value)  / 10 # N -> kg (g=10 m/s^2 dla przyspieszenia ziemskiego)

                res_paras.loc[i] = pd.Series({'sigvm':sigvm,'sfz':sfz})
            else:
                print(f"Obliczenia nie powiodły się dla mutacji {i} z parametrami {[x for x in para_df.loc[i]]}")
                res_paras.loc[i] = pd.Series({'sigvm':np.nan,'sfz':np.nan})

# Połączenie parametrów i wyników do eksportu i wykreślenia
para_out = pd.concat([para_dfo,res_paras], axis=1)

# ---------------------- Eksport wyników ---------------------
para_out.to_excel('./ParaS-VMStress_Thick-min_results.xlsx')

# --------------------- Wykresy ---------------------
fig, ax = plt.subplots(figsize=(12, 9))
ax.scatter(para_out['sfz'], para_out['sigvm'], color='tab:blue', s=50, label='Mutacje')
y_offset = (para_out['sigvm'].max() - para_out['sigvm'].min()) * 0.03
for x, y, mut, t1, t2 in zip(para_out['sfz'], para_out['sigvm'], para_out.index, para_out['t_1'], para_out['t_2']):
    ax.text(x , y+y_offset, f"M{int(mut)}, {t1:.0f}, {t2:.0f} mm", fontsize=8, va='top', ha='center')
ax.axhline(y=f_y, color='red', linestyle='--', linewidth=1, label=f'Granica plastyczności ({f_y} MPa)')
ax.set_title('Naprężenie vs. Masa dla mutacji parametrów')
ax.set_xlabel('Masa [kg]')
ax.set_ylabel('Naprężenie von Misesa [MPa]')
ax.legend()
fig.tight_layout()
fig.savefig('./ParaS-VMStress_Thick-min_results.png', dpi=200)
plt.show()

"""ParaS-LBA_MeshConvergence-min.py

Powiązany model RFEM 6: "MeshSensitivity-LBA_AlC.rf6"
Model: https://www.dlubal.com/en/downloads-and-information/examples-and-tutorials/models-to-download/006080

Przeprowadzenie badania zbieżności siatki w RFEM poprzez zmianę parametru globalnego
`l_fe`, zebranie krytycznego współczynnika obciążenia, porównanie wariantów, zapisanie i wykreślenie wyników.

Sekcje:
- Importy
- Konfiguracja
- Wykonanie główne
- Wykresy
"""

# --------------------- Importy ---------------------
# Biblioteka standardowa
import time

# Strony trzecie
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from dlubal.api import rfem

# ------------------- Konfiguracja -------------------
FE_SIZES_MM = [15, 12, 10, 9, 8, 7, 6, 5, 4, 3] # Lista rozmiarów elementów skończonych do przeszukania w milimetrach
PARAM_NAME = "l_fe"                             # Nazwa parametru globalnego w modelu
CONVERGENCE_THRESHOLD_PCT = 1.0                 # Próg zbieżności w procentach
FCR_ANA = 1065                                  # Analityczny krytyczny współczynnik obciążenia (dla odniesienia na wykresie)

# ------------------- Funkcje -------------------
def set_glpa(dl_app, p_name, p_value,):
    """Ustawia parametr globalny p_name na określoną wartość p_value"""
    params = dl_app.get_object_list([rfem.global_parameters.GlobalParameter()])
    for p in params:
        if p.name == p_name:
            p.value = p_value
            dl_app.update_object(p)
            check=True
            return True
    if not check:
        raise ValueError(f"  Błąd: Nie znaleziono parametru '{p_name}'!")
    
# ------------------- Wykonanie główne -------------------
with rfem.Application() as rf_app:
    # Sprawdź połączenie i wyświetl informacje o modelu
    app_info = rf_app.get_application_info()
    print("Informacje o aplikacji:", app_info)

    results = pd.DataFrame(columns=['l_fe_mm', 'n_elements', 'c_time_s', 'f_cr', 'delta_pct'])
    f_prev = None
    i = 0
    for l_fe in FE_SIZES_MM:
        i += 1
        l_fe_m = l_fe / 1000.0  # Konwersja mm na m do użycia w RFEM
        if not set_glpa(dl_app=rf_app, p_name=PARAM_NAME, p_value=l_fe_m):
            break
        # Ponowna generacja siatki
        rf_app.delete_mesh()
        rf_app.generate_mesh(skip_warnings=True)
        n_elem = rf_app.get_mesh_statistics().surface_2D_finite_elements
        # Uruchom obliczenia
        t0 = time.perf_counter()
        calculation = rf_app.calculate_all(skip_warnings=True)
        t1 = time.perf_counter()
        c_time = t1 - t0
        # Wyniki
        if calculation.succeeded:
            f_cr = rf_app.get_results(
            results_type=rfem.results.STABILITY_ANALYSIS_CRITICAL_LOAD_FACTORS,
            filters=[
            rfem.results.ResultsFilter(column_id='loading', filter_expression='LC1'),
            ]      
            ).data.loc[0].f

            if f_prev is not None and f_prev != 0 and f_cr is not None:
                delta_pct = abs(f_cr - f_prev) / abs(f_prev) * 100.0
                delta_str = f"{delta_pct:>4.2f}"
            else:
                delta_pct = None
                delta_str = f"{'---':>4}"

            if f_cr is not None:
                print(f"Mutacja {i}: l_fe = {l_fe:.1f} mm | n_elements = {n_elem} | f_cr = {f_cr:.2f} | delta = {delta_str} %")
            results.loc[i] = pd.Series({
                'l_fe_mm': l_fe, 'n_elements': n_elem, 'c_time_s': c_time, 
                'f_cr': f_cr, 'delta_pct': delta_pct
                })
            f_prev = f_cr
        else:
            print(f"Obliczenia nie powiodły się dla mutacji {i} z rozmiarem siatki {l_fe:.1f} mm")
            results.loc[i] = pd.Series({
                'l_fe_mm': l_fe, 'n_elements': n_elem, 'c_time_s': np.nan, 
                'f_cr': np.nan, 'delta_pct': np.nan
                })

# Oblicz 1 / n_elements dla drugiego wykresu
results['inv_n_elements'] = 1.0 / results['n_elements']
# Oblicz względne odchylenie od minimalnego krytycznego współczynnika obciążenia
results['rel_dev_f_cr_min'] = (results['f_cr'] - results['f_cr'].min()) / results['f_cr'].min() * 100

# Analiza zbieżności
print(f"\n  Analiza zbieżności (próg: < {CONVERGENCE_THRESHOLD_PCT} % zmiany):")
convergence_point = None
for (l_mm, n, a, d) in results.loc[:, ['l_fe_mm', 'n_elements', 'f_cr', 'delta_pct']].itertuples(index=False):
    if d is not None and d < CONVERGENCE_THRESHOLD_PCT:
        convergence_point = (l_mm, n, a)
        break
if convergence_point:
    print(f"  Zbieżność osiągnięta przy l_fe = {convergence_point[0]:.1f} mm")
    print(f"  f_cr = {convergence_point[2]:.2f}  |  n_elements = {convergence_point[1]}")
else:
    print("  Zbieżność nie została jeszcze osiągnięta – wymagana jest drobniejsza siatka.")

# ---------------------- Eksport wyników ---------------------
results.to_excel('./ParaS-LBA_MeshConvergence-min_results.xlsx')

# ------------------- Wykresy -------------------
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Wykres 1: f_cr vs. n_elements
ax1 = axes[0, 0]
ax1.plot(results['n_elements'], results['f_cr'], 'o-', color='tab:blue', markersize=6)
ax1.axhline(y=FCR_ANA, color='red', linestyle='--', linewidth=1, label=f'$F_{{cr,ana}}$ = {FCR_ANA} kN')
ax1.set_xlabel('Liczba elementów')
ax1.set_ylabel('Krytyczny współczynnik obciążenia [kN]')
ax1.set_title('Krytyczny współczynnik obciążenia vs. Liczba elementów')
ax1.grid(True, alpha=0.3)
ax1.legend()

# Wykres 2: f_cr vs. 1/n_elements
ax2 = axes[0, 1]
ax2.plot(results['inv_n_elements'], results['f_cr'], 'o-', color='tab:blue', markersize=6)
ax2.axhline(y=FCR_ANA, color='red', linestyle='--', linewidth=1, label=f'$F_{{cr,ana}}$ = {FCR_ANA} kN')
ax2.set_xlabel('1 / Liczba elementów')
ax2.set_ylabel('Krytyczny współczynnik obciążenia [kN]')
ax2.set_title('Krytyczny współczynnik obciążenia vs. 1 / Liczba elementów')
ax2.grid(True, alpha=0.3)
ax2.legend()

# Wykres 3: Względne odchylenie vs. Czas obliczeń
ax3 = axes[1, 0]
ax3.plot(results['c_time_s'], results['rel_dev_f_cr_min'], 'o-', color='tab:blue', markersize=6)
ax3.set_xlabel('Czas obliczeń [s]')
ax3.set_ylabel('Względne odchylenie [%]')
ax3.set_title('Względne odchylenie vs. Czas obliczeń')
ax3.grid(True, alpha=0.3)

# Wykres 4: Względne odchylenie vs. l_fe_mm
ax4 = axes[1, 1]
ax4.plot(results['l_fe_mm'], results['rel_dev_f_cr_min'], 'o-', color='tab:blue', markersize=6)
ax4.set_xlabel('Rozmiar siatki ES [mm]')
ax4.set_ylabel('Względne odchylenie [%]')
ax4.set_title('Względne odchylenie vs. Rozmiar siatki ES')
ax4.grid(True, alpha=0.3)

fig.suptitle('Analiza zbieżności siatki', fontsize=14, fontweight='bold')
fig.tight_layout()
fig.savefig('./ParaS-LBA_MeshConvergence-min_results.png', dpi=200)
plt.show()