Skip to content

Commit

Permalink
add to BUG model generation
Browse files Browse the repository at this point in the history
  • Loading branch information
@ACea15
ACea15 committed Oct 8, 2024
1 parent 983b213 commit 3f47e8e
Show file tree
Hide file tree
Showing 3 changed files with 139 additions and 34 deletions.
165 changes: 136 additions & 29 deletions examples/BUG/modelgeneration.org
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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 #
Expand Down Expand Up @@ -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
Expand All @@ -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_<<parameters_modal(output=sol)>>.bdf",
f"BUG103_<<parameters_modal(output=sol)>>.bdf",
scale = 100.)
#+end_src

* DLM generation
Expand Down Expand Up @@ -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]
Expand Down Expand Up @@ -395,7 +416,6 @@

dlm_gafs.build_model()
dlm_gafs.model.write_bdf(f".NASTRAN/GAFs/{gaf_label}.bdf")

#+end_src

** Build Nastran
Expand Down Expand Up @@ -457,13 +477,100 @@
INCLUDE ./DLMs/<<parameters_dlm(output="label_dlm")>>.bdf
INCLUDE ./gafs/<<parameters_gafs(output="label_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<<parameters_gafs(output="label_gaf")>>-<<parameters_gafs(output="num_modes")>>.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<<parameters_gafs(output="label_gaf")>>-<<parameters_gafs(output="num_modes")>>.op4")
op4_Qhj = op4.read_op4(f"./NASTRAN/data_out/Qhh<<parameters_gafs(output="label_gaf")>>-<<parameters_gafs(output="num_modes")>>.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)
2 changes: 1 addition & 1 deletion feniax/aeromodal/roger.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 2 additions & 4 deletions feniax/unastran/op2reader.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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})
Expand All @@ -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())
Expand Down

0 comments on commit 3f47e8e

Please sign in to comment.