diff --git a/examples/BUG/modelgeneration.org b/examples/BUG/modelgeneration.org index 35e245a..dd5d186 100644 --- a/examples/BUG/modelgeneration.org +++ b/examples/BUG/modelgeneration.org @@ -152,8 +152,21 @@ ENDDATA #+end_src +#+begin_src python :results none :tangle generate_aero.py + import feniax.unastran.op2reader as op2reader + # op2 = op2reader.NastranReader(op2name=f"./data_out/Phi{num_modes}") + op2 = OP2() + op2.read_op2("") + eig1 = op2.eigenvectors[1] + eigen = eig1.data + SAVE = True + # if SAVE: + # np.save(f"../FEM/eigenvecs_{num_modes}.npy", modes4sims.T) + # np.save(f"../FEM/eigenvals_{num_modes}.npy", eigs) +#+end_src + ** Free, asets, OP2 (eao) -#+begin_src org :tangle "./NASTRAN/BUG103_cao.bdf" :noweb yes +#+begin_src org :tangle "./NASTRAN/BUG103_eao.bdf" :noweb yes SOL 103 CEND TITLE=BUG model # @@ -198,6 +211,18 @@ ENDDATA #+end_src +#+begin_src python :results none :tangle generate_aero.py + import feniax.unastran.op2reader as op2reader + # op2 = op2reader.NastranReader(op2name=f"./data_out/Phi{num_modes}") + op2 = OP2() + op2.read_op2("") + eig1 = op2.eigenvectors[1] + eigen = eig1.data + SAVE = True + # if SAVE: + # np.save(f"../FEM/eigenvecs_{num_modes}.npy", modes4sims.T) + # np.save(f"../FEM/eigenvals_{num_modes}.npy", eigs) +#+end_src ** Run Nastran #+begin_src bash :results none @@ -207,41 +232,36 @@ #+end_src ** Build modes in OP4 and in ASETs -#+begin_src python :session py1 :results none :tangle generate_aero.py +#+begin_src python :results none :tangle generate_aero.py #from feniax.utils import write_op4modes import feniax.unastran.op4handler as op4handler - eigs, modes = op4handler.write_op4modes("./run_caof", - num_modes, - op4_name=f"./data_out/Phi{num_modes}", - return_modes=True) + import numpy as np + eigs, modes = op4handler.write_op4modes(f"./data_out/Phi{num_modes}", + num_modes, + #op4_name=f"./data_out/Phi{num_modes}", + return_modes=True) bdf = BDF() bdf.read_bdf(self.bdf_file) node_ids = bdf.node_ids - sorted_nodeids = sorted(node_ids) - - # num_nodes = 79 - # eigs = np.array(eigs) - # modes4sims = np.zeros((6*(num_nodes - 1), num_modes)) - # for i in range(num_modes): - # modes4sims[:, i] = np.hstack(modes[i,:(num_nodes - 1)]) - # SAVE = True - # #np.load("../FEM/Ka.npy") - # #scipy.linalg.eigh(Ka, Ma) - # if SAVE: - # np.save("../FEM/eigenvecs.npy", modes4sims) - # np.save("../FEM/eigenvals.npy", eigs) + asets_ids = bdf.asets[0].node_ids + asets_idsfull = np.array([sorted_nodeids.index(ai) for ai in asets_ids]) + modes4simulations = modes[:, asets_idsfull] + SAVE = True + if SAVE: + np.save(f"../FEM/eigenvecs_{num_modes}.npy", modes4sims.T) + np.save(f"../FEM/eigenvals_{num_modes}.npy", eigs) #+end_src - - ** VTK modes #+begin_src python :session py1 :results none :tangle py1.py import feniax.plotools.nastranvtk.modes as modes import importlib importlib.reload(modes) - modes.vtk_fromop2("./SOL103/run_cao.bdf", "./SOL103/run_cao.op2", scale = 100.) + modes.vtk_fromop2(f"BUG103_<>.bdf", + f"BUG103_<>.bdf", + scale = 100.) #+end_src * DLM generation @@ -355,7 +375,8 @@ import numpy as np import feniax.unastran.aero as nasaero import pickle - + sol = "cao" + num_modes = 50 mach = 0.8 Mach = str(mach).replace('.','_') machs = [mach] @@ -395,7 +416,6 @@ dlm_gafs.build_model() dlm_gafs.model.write_bdf(f".NASTRAN/GAFs/{gaf_label}.bdf") - #+end_src ** Build Nastran @@ -457,13 +477,100 @@ INCLUDE ./DLMs/<>.bdf INCLUDE ./gafs/<>.bdf #+end_src -** Read GAFs -#+begin_src python :results none :noweb yes :tangle read_gafs.py +** Roger RFA +#+begin_src python :results none :noweb yes :tangle roger.py import pyNastran.op4.op4 as op4 - Qhh = op4.read_op4(f"Qhh<>-<>.op4") + import feniax.aeromodal.roger as roger + import plotly.express as px + import plotly.graph_objects as go + + op4_Qhh = op4.read_op4(f"./NASTRAN/data_out/Qhh<>-<>.op4") + op4_Qhj = op4.read_op4(f"./NASTRAN/data_out/Qhh<>-<>.op4") + + try: + qhh = jnp.array(op4_Qhh["Q_HH"].data) + qhj = jnp.array(op4_Qhj["Q_HJ"].data) + except AttributeError: + qhh = jnp.array(op4_Qhh["Q_HH"][1]) + qhj = jnp.array(op4_Qhj["Q_HJ"][1]) + + poles = jnp.linspace(1e-3, 1, 40) + qhj_new = roger.stackQk_realimag(qhj) + qhh_new = roger.stackQk_realimag(qhh) + + k_matrix = roger.frequency_matrix(reduced_freqs[1:], poles) + roger_matricesQhj = roger.rogerRFA(k_matrix, qhj_new) + roger_matricesQhh = roger.rogerRFA(k_matrix, qhh_new) + + ks = jnp.hstack([1e-6, k_array]) + Qk_hj = roger.Q_RFA(reduced_freqs, roger_matricesQhj, poles) + Qk_hh = roger.Q_RFA(reduced_freqs, roger_matricesQhh, poles) + + PLOT = True + if PLOT: + + i = 1 + j = 50 + #i = 20 + #j = 18 + + fig = go.Figure() + # fig.add_trace(go.Scatter(x=qhh[:,i,j].real, y=qhh[:,i,j].imag), + # mode='makers', + # # name='lines' + # ) + fig.add_trace( + go.Scatter( + x=qhh[:, i, j].real, + y=qhh[:, i, j].imag, + mode="markers", + # name='lines' + ), + ) + + fig.add_trace( + go.Scatter( + x=Qk_hh[:, i, j].real, + y=Qk_hh[:, i, j].imag, + mode="lines", + # name='lines' + ), + ) + # fig.add_trace(px.scatter(x=qhh[:,i,j].real, y=qhh[:,i,j].imag)) + + # fig = px.scatter(x=qhh[:,i,j].real, y=qhh[:,i,j].imag) + fig.show() + + # Q_RFA(ki, roger_matrices, poles) + + fig = go.Figure() + # fig.add_trace(go.Scatter(x=qhh[:,i,j].real, y=qhh[:,i,j].imag), + # mode='makers', + # # name='lines' + # ) + fig.add_trace( + go.Scatter( + x=qhj[:, i, j].real, + y=qhj[:, i, j].imag, + mode="markers", + # name='lines' + ), + ) + + fig.add_trace( + go.Scatter( + x=Qk_hj[:, i, j].real, + y=Qk_hj[:, i, j].imag, + mode="lines", + # name='lines' + ), + ) + # fig.add_trace(px.scatter(x=qhh[:,i,j].real, y=qhh[:,i,j].imag)) + + # fig = px.scatter(x=qhh[:,i,j].real, y=qhh[:,i,j].imag) + fig.show() + #+end_src -** Roger RFA - * Gust solution (146) diff --git a/feniax/aeromodal/roger.py b/feniax/aeromodal/roger.py index e077f24..f7d29da 100644 --- a/feniax/aeromodal/roger.py +++ b/feniax/aeromodal/roger.py @@ -166,9 +166,9 @@ def Q_RFAki(ki, roger_matrices, poles): return Qk +Q_RFA = jax.vmap(Q_RFAki, in_axes=(0, None, None)) if __name__ == "__main__": - Q_RFA = jax.vmap(Q_RFAki, in_axes=(0, None, None)) def err_ki(ki, Qki_dlm, roger_matrices, poles, order=None): Qki_roger = Q_RFA(ki, roger_matrices, poles) diff --git a/feniax/unastran/op2reader.py b/feniax/unastran/op2reader.py index 66f3b6f..61ee1ec 100644 --- a/feniax/unastran/op2reader.py +++ b/feniax/unastran/op2reader.py @@ -13,8 +13,8 @@ def __init__(self, op2name=None, bdfname=None, static=None): self.op2name = op2name self.bdfname = bdfname self.static = static - - def readModel(self, asets_only=False): + self._read + def _readModel(self, asets_only=False): if self.op2name: self.op2 = OP2() self.op2.set_additional_matrices_to_read({b"OPHP": False, b"OUG1": False}) @@ -36,8 +36,6 @@ def eigenvectors(self): eig1 = self.op2.eigenvectors[1] eigen = eig1.data NumNodes = np.shape(eig1.data[0])[0] - # model.eigenvectors[1].modes - # model.eigenvectors[1]._times def displacements(self): subcases = sorted(self.op2.displacements.keys())