Pynastranを用いたOP2ファイルからの結果出力

LLMで情報収集中

Pynastranで結果の取得をしたい。用件を箇条書きで示します。▪︎計算はノーマルモード解析▪︎プロパティに含まれる要素ごとに合計値を出力▪︎出力したい結果はひずみエネルギー▪︎最終的にはモデル全体のひずみエネルギーの合計値とプロパティ単位の合計値の比を求め、プロパティ単位での全体ひずみエネルギーに対する割合(寄与度)を取得したい。必要情報があれば先に聞いてください。ネット検索必要なら利用してください。

了解です。pyNastranで SOL103（ノーマルモード）の要素ひずみエネルギー（ESE）を読み、Property（PID）単位に合計 → 全体合計に対する寄与度（比率）を出す、という流れで組めます。 --- まず前提：NASTRAN側で「ESE」を出力しておく OP2やPCHに ESE が出ていないと、pyNastran側では取得できません。一般に Case Control に例えば下記のような指定を入れます（GUIだと「Strain Energy」をONにするのと同義）。 ESE(PUNCH)=ALL（PCH出力）（OP2で欲しいなら）解析のI/O設定で OP2 を有効化（環境によって PARAM,POST,-1 等が必要な場合あり） ※ OP2内では「モードの要素ひずみエネルギー」系テーブル（例：ONRGY1/ONRGY2）として扱われます。 --- pyNastranでやること（設計） 1. BDFを読む：eid -> pid の対応表を作る（要素からPIDを取る） 2. OP2を読む：strain_energy（ESE）結果を取り出すpyNastranは op2.op2_results.strain_energy.cquad4_strain_energy など「要素タイプ別」に保持します。 3. モードごとに要素ひずみエネルギーを eid で取り出し pid に集約して合計全体合計で割って寄与度（割合）算出 4. CSV等に出力 --- サンプルコード（OP2 + BDFから、PID別寄与度を出す） > なるべく「要素タイプが混在しても動く」ように、strain_energy配下を introspection して総なめします。また、ESEの列は「energy / percent / density」のような構成になり得るので（古いWikiにもその記述があります）、基本は get_headers() が取れればそれに従い、取れない場合は先頭列を energy とみなす、で堅めにしています。 from __future__ import annotations from collections import defaultdictfrom pathlib import Path import numpy as npimport pandas as pd from pyNastran.bdf.bdf import BDFfrom pyNastran.op2.op2 import OP2 def build_eid_to_pid_map(bdf: BDF) -> dict[int, int]: """要素EID -> Property ID(PID) の対応を作る。PIDを持たない要素は除外。""" eid_to_pid: dict[int, int] = {} for eid, elem in bdf.elements.items(): # RBE2等 pid を持たない要素もある pid = getattr(elem, "pid", None) if pid is None: continue try: pid_int = int(pid) except Exception: continue eid_to_pid[int(eid)] = pid_int return eid_to_pid def iter_strain_energy_result_objects(op2: OP2): """ op2.op2_results.strain_energy 配下のテーブル(dict)を列挙する。例: cquad4_strain_energy, ctetra_strain_energy ... など。4 """ se = op2.op2_results.strain_energy for name in dir(se): if not name.endswith("_strain_energy"): continue obj = getattr(se, name, None) if isinstance(obj, dict) and obj: # subcase -> result_object のdict yield name, obj def pick_energy_column(result_obj) -> int: """ result_obj.data のどの列が strain energy かを推定する。 headersが取れるなら、それを優先。取れないなら 0 列目を採用。 """ if hasattr(result_obj, "get_headers"): headers = result_obj.get_headers() # 例: ['strain_energy', 'percent', 'strain_energy_density'] のようなケースを想定 for i, h in enumerate(headers): h_low = str(h).lower() if "strain" in h_low and "energy" in h_low and "percent" not in h_low and "density" not in h_low: return i if h_low in ("energy", "strain_energy"): return i # それでも見つからなければ、percent/density以外の最初の列を採用 for i, h in enumerate(headers): h_low = str(h).lower() if "percent" in h_low or "density" in h_low: continue return i return 0 def compute_property_contribution( bdf_path: str | Path, op2_path: str | Path, subcase: int | None = None, modes: list[int] | None = None, # 0-based index（pyNastranのtime/mode index想定）) -> pd.DataFrame: # --- read BDF (for eid->pid) --- bdf = BDF(debug=False) bdf.read_bdf(str(bdf_path), xref=False) eid_to_pid = build_eid_to_pid_map(bdf) # --- read OP2 (for strain energy) --- op2 = OP2(debug=False) op2.read_op2(str(op2_path)) rows = [] for table_name, subcase_dict in iter_strain_energy_result_objects(op2): for key, result in subcase_dict.items(): # subcase指定がある場合：keyが int とは限らない（タプルキーになる例あり）5 # ここは「intなら一致」「タプルなら先頭のsubcaseを使う」方針で吸収 if subcase is not None: k_subcase = None if isinstance(key, int): k_subcase = key elif isinstance(key, tuple) and len(key) >= 1 and isinstance(key[0], int): k_subcase = key[0] if k_subcase != subcase: continue # result.element : (nelements,) # result.data : (nmodes, nelements, ncols) を想定 eids = np.asarray(getattr(result, "element", []), dtype=int) if eids.size == 0 or not hasattr(result, "data"): continue data = np.asarray(result.data) if data.ndim != 3: continue energy_col = pick_energy_column(result) nmodes = data.shape[0] mode_indices = range(nmodes) if modes is None else [m for m in modes if 0 <= m < nmodes] for mi in mode_indices: energy = data[mi, :, energy_col].astype(float) pid_sum = defaultdict(float) total = 0.0 for eid, e in zip(eids, energy): pid = eid_to_pid.get(int(eid)) if pid is None: continue pid_sum[pid] += float(e) total += float(e) if total == 0.0: continue for pid, se_val in pid_sum.items(): rows.append( { "table": table_name, "subcase_key": str(key), "mode_index": mi, "pid": pid, "strain_energy": se_val, "contribution": se_val / total, "total_strain_energy_mode": total, } ) df = pd.DataFrame(rows) if df.empty: return df # 見やすく：寄与度[%] df["contribution_pct"] = 100.0 * df["contribution"] # 典型的には table が要素タイプごとに分かれてるので、 # 「PID単位に最終集約（要素タイプを跨いで合算）」した集計表も出しておく agg = ( df.groupby(["subcase_key", "mode_index", "pid"], as_index=False) .agg( strain_energy=("strain_energy", "sum"), total_strain_energy_mode=("total_strain_energy_mode", "max"), ) ) agg["contribution"] = agg["strain_energy"] / agg["total_strain_energy_mode"] agg["contribution_pct"] = 100.0 * agg["contribution"] agg = agg.sort_values(["subcase_key", "mode_index", "contribution_pct"], ascending=[True, True, False]) return agg if __name__ == "__main__": # 入力 bdf_path = "model.bdf" op2_path = "model.op2" # 例：subcase=1 の全モードで実行（modesは 0-based index） df = compute_property_contribution(bdf_path, op2_path, subcase=1, modes=None) print(df.head(30)) df.to_csv("strain_energy_contribution_by_pid.csv", index=False) --- ここまでで満たせる要件 ✅ ノーマルモード解析（SOL103）のESE結果を取得（OP2に入っている前提） ✅ PID（プロパティ）単位で合計値 ✅ モデル全体の合計に対する PIDごとの比（寄与度） ✅ CSV出力まで --- 先に確認したい「必要情報」（最小限）あなたの環境に合わせてコードの“ツメ”を調整したいので、次だけ教えてください（1行でOKです）： 1. OP2は出ていますか？（ESEを入れて再計算済みか） 2. 使用ソルバは MSC / NX / OptiStruct のどれですか？（OP2の癖対策用。ESEでキーがタプル化する例もあります） 3. 寄与度は「各モードごと」が欲しいですか？それとも「指定モード範囲で総和」（例：1〜6次モード合算）ですか？この3点が分かれば、あなたの出したい最終形式（例：モード別シート、合算1枚、PID名も併記…など）に合わせて完成形に寄せます。

OP2にひずみエネルギーは出力済み。ソルバはNX。寄与度は各モードごとで欲しいです。

了解です（NX／OP2にESE出力済み／寄与度は各モードごと）。 NX Nastran の ESE（要素ひずみエネルギー）を pyNastran で読むときの “ハマりどころ” は主にこれです： strain_energy の辞書キーが int（=subcase番号）ではなく、tuple になることがある（例：(1, 2, 1, 0, 0, '', '') みたいなキー）→ 「subcase=1だけ欲しい」のに単純に dict[subcase] ができないことがあるので、キーから subcase を抽出してフィルタします。以下に、NXのtupleキーも吸収して、(モード×PID) の寄与度をCSVにするコードを完成形として置きます。

from __future__ import annotations

from collections import defaultdict
from pathlib import Path
import numpy as np
import pandas as pd

from pyNastran.bdf.bdf import BDF
from pyNastran.op2.op2 import OP2


def build_eid_to_pid_map(bdf: BDF) -> dict[int, int]:
    """要素EID -> PID。pidを持たない要素(RBE等)は除外。"""
    out: dict[int, int] = {}
    for eid, elem in bdf.elements.items():
        pid = getattr(elem, "pid", None)
        if pid is None:
            continue
        try:
            out[int(eid)] = int(pid)
        except Exception:
            continue
    return out


def get_subcase_from_key(key) -> int | None:
    """
    NXでありがちな strain_energy のキー（tuple）から subcase を取り出す。
    例: (1, 2, 1, 0, 0, '', '') の先頭が isubcase。1
    """
    if isinstance(key, int):
        return key
    if isinstance(key, tuple) and len(key) >= 1 and isinstance(key[0], int):
        return key[0]
    return None


def iter_strain_energy_tables(op2: OP2):
    """
    op2.op2_results.strain_energy 配下の *_strain_energy dict を列挙。
    要素タイプ別（cquad4/ctetra等）に分かれている。2
    """
    se = op2.op2_results.strain_energy
    for name in dir(se):
        if not name.endswith("_strain_energy"):
            continue
        dct = getattr(se, name, None)
        if isinstance(dct, dict) and dct:
            yield name, dct


def pick_energy_column(result_obj) -> int:
    """dataのどの列が strain energy か推定（headersがあれば優先）。"""
    if hasattr(result_obj, "get_headers"):
        headers = [str(h).lower() for h in result_obj.get_headers()]
        # まず "strain" & "energy" を含み、percent/density ではない列を探す
        for i, h in enumerate(headers):
            if ("strain" in h and "energy" in h) and ("percent" not in h) and ("density" not in h):
                return i
        # 次点：energy/strain_energy
        for i, h in enumerate(headers):
            if h in ("energy", "strain_energy"):
                return i
        # それでもダメなら percent/density以外の最初
        for i, h in enumerate(headers):
            if ("percent" in h) or ("density" in h):
                continue
            return i
    return 0


def get_mode_numbers(result_obj, nmodes: int) -> list[int]:
    """
    pyNastran側の「モード番号（1始まり）」を可能なら取得し、
    取れなければ 1..nmodes を採用。
    """
    # 多くの結果オブジェクトは _times に固有値/周波数/モード番号が入ることがある
    modes = getattr(result_obj, "modes", None)
    if modes is not None and len(modes) == nmodes:
        try:
            return [int(x) for x in modes]
        except Exception:
            pass

    times = getattr(result_obj, "_times", None)
    if times is not None and len(times) == nmodes:
        # SOL103だと "mode number" が入る場合と固有値が入る場合があるので、
        # ここは「整数っぽければそれを採用」くらいに留める
        try:
            as_int = [int(x) for x in times]
            # 1..nmodesに近いならモード番号扱い
            if set(as_int) == set(range(1, nmodes + 1)):
                return as_int
        except Exception:
            pass

    return list(range(1, nmodes + 1))


def compute_modal_pid_contribution(
    bdf_path: str | Path,
    op2_path: str | Path,
    subcase: int = 1,
) -> pd.DataFrame:
    # --- read BDF ---
    bdf = BDF(debug=False)
    bdf.read_bdf(str(bdf_path), xref=False)
    eid_to_pid = build_eid_to_pid_map(bdf)

    # --- read OP2 ---
    op2 = OP2(debug=False)
    op2.read_op2(str(op2_path))

    # --- accumulate (mode, pid) across all element types ---
    # NXでは element-type別テーブルに分かれるので、最後に跨って合算する
    accum = defaultdict(float)  # (mode_no, pid) -> sum_energy
    totals = defaultdict(float) # mode_no -> total_energy

    for table_name, dct in iter_strain_energy_tables(op2):
        for key, res in dct.items():
            if get_subcase_from_key(key) != subcase:
                continue

            eids = np.asarray(getattr(res, "element", []), dtype=int)
            if eids.size == 0 or not hasattr(res, "data"):
                continue

            data = np.asarray(res.data)
            if data.ndim != 3:
                continue

            nmodes = data.shape[0]
            energy_col = pick_energy_column(res)
            mode_nos = get_mode_numbers(res, nmodes)

            # data: (nmodes, nelems, ncols)
            for mi, mode_no in enumerate(mode_nos):
                energy = data[mi, :, energy_col].astype(float)

                # PID集約
                mode_total = 0.0
                for eid, e in zip(eids, energy):
                    pid = eid_to_pid.get(int(eid))
                    if pid is None:
                        continue
                    accum[(mode_no, pid)] += float(e)
                    mode_total += float(e)

                totals[mode_no] += mode_total

    # --- build dataframe ---
    rows = []
    for (mode_no, pid), se in accum.items():
        total = totals.get(mode_no, 0.0)
        if total == 0.0:
            continue
        rows.append({
            "subcase": subcase,
            "mode_no": mode_no,          # 1始まりのモード番号（推定）
            "pid": pid,
            "strain_energy": se,
            "total_strain_energy": total,
            "contribution": se / total,
            "contribution_pct": 100.0 * (se / total),
        })

    df = pd.DataFrame(rows)
    if df.empty:
        return df

    df = df.sort_values(["mode_no", "contribution_pct"], ascending=[True, False]).reset_index(drop=True)
    return df


if __name__ == "__main__":
    bdf_path = "model.bdf"
    op2_path = "model.op2"
    subcase = 1

    df = compute_modal_pid_contribution(bdf_path, op2_path, subcase=subcase)
    print(df.head(30))
    df.to_csv("modal_strain_energy_contribution_by_pid.csv", index=False)

--- NX向け：SOL103 ESE → PID合計 → 寄与度（各モード） from __future__ import annotations from collections import defaultdictfrom pathlib import Pathimport numpy as npimport pandas as pd from pyNastran.bdf.bdf import BDFfrom pyNastran.op2.op2 import OP2 def build_eid_to_pid_map(bdf: BDF) -> dict[int, int]: """要素EID -> PID。pidを持たない要素(RBE等)は除外。""" out: dict[int, int] = {} for eid, elem in bdf.elements.items(): pid = getattr(elem, "pid", None) if pid is None: continue try: out[int(eid)] = int(pid) except Exception: continue return out def get_subcase_from_key(key) -> int | None: """ NXでありがちな strain_energy のキー（tuple）から subcase を取り出す。例: (1, 2, 1, 0, 0, '', '') の先頭が isubcase。1 """ if isinstance(key, int): return key if isinstance(key, tuple) and len(key) >= 1 and isinstance(key[0], int): return key[0] return None def iter_strain_energy_tables(op2: OP2): """ op2.op2_results.strain_energy 配下の *_strain_energy dict を列挙。要素タイプ別（cquad4/ctetra等）に分かれている。2 """ se = op2.op2_results.strain_energy for name in dir(se): if not name.endswith("_strain_energy"): continue dct = getattr(se, name, None) if isinstance(dct, dict) and dct: yield name, dct def pick_energy_column(result_obj) -> int: """dataのどの列が strain energy か推定（headersがあれば優先）。""" if hasattr(result_obj, "get_headers"): headers = [str(h).lower() for h in result_obj.get_headers()] # まず "strain" & "energy" を含み、percent/density ではない列を探す for i, h in enumerate(headers): if ("strain" in h and "energy" in h) and ("percent" not in h) and ("density" not in h): return i # 次点：energy/strain_energy for i, h in enumerate(headers): if h in ("energy", "strain_energy"): return i # それでもダメなら percent/density以外の最初 for i, h in enumerate(headers): if ("percent" in h) or ("density" in h): continue return i return 0 def get_mode_numbers(result_obj, nmodes: int) -> list[int]: """ pyNastran側の「モード番号（1始まり）」を可能なら取得し、取れなければ 1..nmodes を採用。 """ # 多くの結果オブジェクトは _times に固有値/周波数/モード番号が入ることがある modes = getattr(result_obj, "modes", None) if modes is not None and len(modes) == nmodes: try: return [int(x) for x in modes] except Exception: pass times = getattr(result_obj, "_times", None) if times is not None and len(times) == nmodes: # SOL103だと "mode number" が入る場合と固有値が入る場合があるので、 # ここは「整数っぽければそれを採用」くらいに留める try: as_int = [int(x) for x in times] # 1..nmodesに近いならモード番号扱い if set(as_int) == set(range(1, nmodes + 1)): return as_int except Exception: pass return list(range(1, nmodes + 1)) def compute_modal_pid_contribution( bdf_path: str | Path, op2_path: str | Path, subcase: int = 1,) -> pd.DataFrame: # --- read BDF --- bdf = BDF(debug=False) bdf.read_bdf(str(bdf_path), xref=False) eid_to_pid = build_eid_to_pid_map(bdf) # --- read OP2 --- op2 = OP2(debug=False) op2.read_op2(str(op2_path)) # --- accumulate (mode, pid) across all element types --- # NXでは element-type別テーブルに分かれるので、最後に跨って合算する accum = defaultdict(float) # (mode_no, pid) -> sum_energy totals = defaultdict(float) # mode_no -> total_energy for table_name, dct in iter_strain_energy_tables(op2): for key, res in dct.items(): if get_subcase_from_key(key) != subcase: continue eids = np.asarray(getattr(res, "element", []), dtype=int) if eids.size == 0 or not hasattr(res, "data"): continue data = np.asarray(res.data) if data.ndim != 3: continue nmodes = data.shape[0] energy_col = pick_energy_column(res) mode_nos = get_mode_numbers(res, nmodes) # data: (nmodes, nelems, ncols) for mi, mode_no in enumerate(mode_nos): energy = data[mi, :, energy_col].astype(float) # PID集約 mode_total = 0.0 for eid, e in zip(eids, energy): pid = eid_to_pid.get(int(eid)) if pid is None: continue accum[(mode_no, pid)] += float(e) mode_total += float(e) totals[mode_no] += mode_total # --- build dataframe --- rows = [] for (mode_no, pid), se in accum.items(): total = totals.get(mode_no, 0.0) if total == 0.0: continue rows.append({ "subcase": subcase, "mode_no": mode_no, # 1始まりのモード番号（推定） "pid": pid, "strain_energy": se, "total_strain_energy": total, "contribution": se / total, 出力されるもの CSVに、モードごとに total_strain_energy（モデル全体合計） pidごとの strain_energy 合計 contribution（0〜1） contribution_pct（%）が入ります。 --- NXでよくある注意（重要） tupleキー問題：上の get_subcase_from_key() で吸収しています（NXのESEでよく起きる挙動として報告あり）。「モード番号」表示：pyNastranの結果オブジェクトによっては modes や _times に “モード番号そのもの” が入らず、固有値など別の値が入る場合があります。その場合このコードは安全側で 1..nmodes を mode_no として採番します（寄与度計算自体は問題なし）。 --- 次に、精度を上げるために1点だけ確認（質問というより“チェック項目”）もし「mode_no が 1,2,3… と綺麗に出てない」場合は、あなたのOP2内でモード番号/固有値がどの配列に入っているかを合わせ込みます。その時は、あなたの環境で次を1回だけ実行して出力を貼ってください（10行程度でOK）： from pyNastran.op2.op2 import OP2op2 = OP2()op2.read_op2("model.op2") se = op2.op2_results.strain_energy# 何か1個（例: cquad4_strain_energy）を選んで最初のresultを覗くdct = se.cquad4_strain_energykey0 = list(dct.keys())[0]res = dct[key0]print("key0 =", key0)print("get_headers =", res.get_headers() if hasattr(res,"get_headers") else None)print("has modes =", hasattr(res,"modes"))print("modes =", getattr(res,"modes", None))print("times =", getattr(res,"_times", None))print("data shape =", res.data.shape) （この情報があると、mode_no を確実に OP2 のモード番号に一致させられます。）ここまでで、要件どおり NX SOL103の“各モードごとのPID寄与度” を取れるはずです。