diff --git a/examples/RafaFlyingBeam/FEM/Ka_m1.npy b/examples/RafaFlyingBeam/FEM/Ka_m1.npy new file mode 100644 index 0000000..0e62a42 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/Ka_m1.npy differ diff --git a/examples/RafaFlyingBeam/FEM/Ka_m2.npy b/examples/RafaFlyingBeam/FEM/Ka_m2.npy new file mode 100644 index 0000000..0e62a42 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/Ka_m2.npy differ diff --git a/examples/RafaFlyingBeam/FEM/Ka_m3.npy b/examples/RafaFlyingBeam/FEM/Ka_m3.npy new file mode 100644 index 0000000..44b0c24 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/Ka_m3.npy differ diff --git a/examples/RafaFlyingBeam/FEM/Ka_m4.npy b/examples/RafaFlyingBeam/FEM/Ka_m4.npy new file mode 100644 index 0000000..44b0c24 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/Ka_m4.npy differ diff --git a/examples/RafaFlyingBeam/FEM/Ma_m1.npy b/examples/RafaFlyingBeam/FEM/Ma_m1.npy new file mode 100644 index 0000000..26c9033 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/Ma_m1.npy differ diff --git a/examples/RafaFlyingBeam/FEM/Ma_m2.npy b/examples/RafaFlyingBeam/FEM/Ma_m2.npy new file mode 100644 index 0000000..7202052 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/Ma_m2.npy differ diff --git a/examples/RafaFlyingBeam/FEM/Ma_m3.npy b/examples/RafaFlyingBeam/FEM/Ma_m3.npy new file mode 100644 index 0000000..26c9033 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/Ma_m3.npy differ diff --git a/examples/RafaFlyingBeam/FEM/Ma_m4.npy b/examples/RafaFlyingBeam/FEM/Ma_m4.npy new file mode 100644 index 0000000..7202052 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/Ma_m4.npy differ diff --git a/examples/RafaFlyingBeam/FEM/eigenvals_m1.npy b/examples/RafaFlyingBeam/FEM/eigenvals_m1.npy new file mode 100644 index 0000000..837db96 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/eigenvals_m1.npy differ diff --git a/examples/RafaFlyingBeam/FEM/eigenvals_m2.npy b/examples/RafaFlyingBeam/FEM/eigenvals_m2.npy new file mode 100644 index 0000000..5abb83d Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/eigenvals_m2.npy differ diff --git a/examples/RafaFlyingBeam/FEM/eigenvals_m3.npy b/examples/RafaFlyingBeam/FEM/eigenvals_m3.npy new file mode 100644 index 0000000..c78cfa8 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/eigenvals_m3.npy differ diff --git a/examples/RafaFlyingBeam/FEM/eigenvals_m4.npy b/examples/RafaFlyingBeam/FEM/eigenvals_m4.npy new file mode 100644 index 0000000..c70cf8b Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/eigenvals_m4.npy differ diff --git a/examples/RafaFlyingBeam/FEM/eigenvecs_m1.npy b/examples/RafaFlyingBeam/FEM/eigenvecs_m1.npy new file mode 100644 index 0000000..a0c0da2 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/eigenvecs_m1.npy differ diff --git a/examples/RafaFlyingBeam/FEM/eigenvecs_m2.npy b/examples/RafaFlyingBeam/FEM/eigenvecs_m2.npy new file mode 100644 index 0000000..9d03784 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/eigenvecs_m2.npy differ diff --git a/examples/RafaFlyingBeam/FEM/eigenvecs_m3.npy b/examples/RafaFlyingBeam/FEM/eigenvecs_m3.npy new file mode 100644 index 0000000..4e03c88 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/eigenvecs_m3.npy differ diff --git a/examples/RafaFlyingBeam/FEM/eigenvecs_m4.npy b/examples/RafaFlyingBeam/FEM/eigenvecs_m4.npy new file mode 100644 index 0000000..3f49103 Binary files /dev/null and b/examples/RafaFlyingBeam/FEM/eigenvecs_m4.npy differ diff --git a/examples/RafaFlyingBeam/FEM/structuralGrid_m1 b/examples/RafaFlyingBeam/FEM/structuralGrid_m1 new file mode 100644 index 0000000..3a7161d --- /dev/null +++ b/examples/RafaFlyingBeam/FEM/structuralGrid_m1 @@ -0,0 +1,4 @@ +# VARIABLES = x y z fe_order components +0.0 0.0 0.0 0 rbeam +1.0 0.0 0.0 1 rbeam +-1.0 0.0 0.0 2 lbeam diff --git a/examples/RafaFlyingBeam/FEM/structuralGrid_m2 b/examples/RafaFlyingBeam/FEM/structuralGrid_m2 new file mode 100644 index 0000000..3a7161d --- /dev/null +++ b/examples/RafaFlyingBeam/FEM/structuralGrid_m2 @@ -0,0 +1,4 @@ +# VARIABLES = x y z fe_order components +0.0 0.0 0.0 0 rbeam +1.0 0.0 0.0 1 rbeam +-1.0 0.0 0.0 2 lbeam diff --git a/examples/RafaFlyingBeam/FEM/structuralGrid_m3 b/examples/RafaFlyingBeam/FEM/structuralGrid_m3 new file mode 100644 index 0000000..a4e1d87 --- /dev/null +++ b/examples/RafaFlyingBeam/FEM/structuralGrid_m3 @@ -0,0 +1,4 @@ +# VARIABLES = x y z fe_order components +0.0 0.0 0.0 0 rbeam +0.8660254 0.0 0.5 1 rbeam +-0.866025 0.0 -0.5 2 lbeam diff --git a/examples/RafaFlyingBeam/FEM/structuralGrid_m4 b/examples/RafaFlyingBeam/FEM/structuralGrid_m4 new file mode 100644 index 0000000..a4e1d87 --- /dev/null +++ b/examples/RafaFlyingBeam/FEM/structuralGrid_m4 @@ -0,0 +1,4 @@ +# VARIABLES = x y z fe_order components +0.0 0.0 0.0 0 rbeam +0.8660254 0.0 0.5 1 rbeam +-0.866025 0.0 -0.5 2 lbeam diff --git a/examples/RafaFlyingBeam/NASTRAN/Model1_103op2.bdf b/examples/RafaFlyingBeam/NASTRAN/Model1_103op2.bdf new file mode 100644 index 0000000..7a34661 --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/Model1_103op2.bdf @@ -0,0 +1,11 @@ +SOL 103 +CEND +TITLE=Rafa's flying beam model # +ECHO=NONE +DISPLACEMENT=ALL +METHOD = 900 +BEGIN BULK +EIGRL,900,,,18 +PARAM,POST,-1 +INCLUDE 'model1.bdf' +ENDDATA diff --git a/examples/RafaFlyingBeam/NASTRAN/Model1_103pch.bdf b/examples/RafaFlyingBeam/NASTRAN/Model1_103pch.bdf new file mode 100644 index 0000000..45f200f --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/Model1_103pch.bdf @@ -0,0 +1,11 @@ +SOL 103 +CEND +TITLE=Rafa's flying beam model # +ECHO=NONE +DISPLACEMENT=ALL +METHOD = 900 +BEGIN BULK +EIGRL,900,,,18 +PARAM,EXTOUT,DMIGPCH +INCLUDE 'model1.bdf' +ENDDATA diff --git a/examples/RafaFlyingBeam/NASTRAN/Model2_103op2.bdf b/examples/RafaFlyingBeam/NASTRAN/Model2_103op2.bdf new file mode 100644 index 0000000..49f3121 --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/Model2_103op2.bdf @@ -0,0 +1,11 @@ +SOL 103 +CEND +TITLE=Rafa's flying beam model # +ECHO=NONE +DISPLACEMENT=ALL +METHOD = 900 +BEGIN BULK +EIGRL,900,,,18 +PARAM,POST,-1 +INCLUDE 'model2.bdf' +ENDDATA diff --git a/examples/RafaFlyingBeam/NASTRAN/Model2_103pch.bdf b/examples/RafaFlyingBeam/NASTRAN/Model2_103pch.bdf new file mode 100644 index 0000000..0ea04db --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/Model2_103pch.bdf @@ -0,0 +1,11 @@ +SOL 103 +CEND +TITLE=Rafa's flying beam model # +ECHO=NONE +DISPLACEMENT=ALL +METHOD = 900 +BEGIN BULK +EIGRL,900,,,18 +PARAM,EXTOUT,DMIGPCH +INCLUDE 'model2.bdf' +ENDDATA diff --git a/examples/RafaFlyingBeam/NASTRAN/Model3_103op2.bdf b/examples/RafaFlyingBeam/NASTRAN/Model3_103op2.bdf new file mode 100644 index 0000000..19ec24b --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/Model3_103op2.bdf @@ -0,0 +1,11 @@ +SOL 103 +CEND +TITLE=Rafa's flying beam model # +ECHO=NONE +DISPLACEMENT=ALL +METHOD = 900 +BEGIN BULK +EIGRL,900,,,18 +PARAM,POST,-1 +INCLUDE 'model3.bdf' +ENDDATA diff --git a/examples/RafaFlyingBeam/NASTRAN/Model3_103pch.bdf b/examples/RafaFlyingBeam/NASTRAN/Model3_103pch.bdf new file mode 100644 index 0000000..f520cd3 --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/Model3_103pch.bdf @@ -0,0 +1,11 @@ +SOL 103 +CEND +TITLE=Rafa's flying beam model # +ECHO=NONE +DISPLACEMENT=ALL +METHOD = 900 +BEGIN BULK +EIGRL,900,,,18 +PARAM,EXTOUT,DMIGPCH +INCLUDE 'model3.bdf' +ENDDATA diff --git a/examples/RafaFlyingBeam/NASTRAN/Model4_103op2.bdf b/examples/RafaFlyingBeam/NASTRAN/Model4_103op2.bdf new file mode 100644 index 0000000..114bf59 --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/Model4_103op2.bdf @@ -0,0 +1,11 @@ +SOL 103 +CEND +TITLE=Rafa's flying beam model # +ECHO=NONE +DISPLACEMENT=ALL +METHOD = 900 +BEGIN BULK +EIGRL,900,,,18 +PARAM,POST,-1 +INCLUDE 'model4.bdf' +ENDDATA diff --git a/examples/RafaFlyingBeam/NASTRAN/Model4_103pch.bdf b/examples/RafaFlyingBeam/NASTRAN/Model4_103pch.bdf new file mode 100644 index 0000000..3c56192 --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/Model4_103pch.bdf @@ -0,0 +1,11 @@ +SOL 103 +CEND +TITLE=Rafa's flying beam model # +ECHO=NONE +DISPLACEMENT=ALL +METHOD = 900 +BEGIN BULK +EIGRL,900,,,18 +PARAM,EXTOUT,DMIGPCH +INCLUDE 'model4.bdf' +ENDDATA diff --git a/examples/RafaFlyingBeam/NASTRAN/model1.bdf b/examples/RafaFlyingBeam/NASTRAN/model1.bdf new file mode 100644 index 0000000..d219edb --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/model1.bdf @@ -0,0 +1,23 @@ +$pyNastran: version=msc +$pyNastran: punch=True +$pyNastran: encoding=utf-8 +$pyNastran: nnodes=3 +$pyNastran: nelements=2 +$NODES +GRID 1 0. 0. 0. +GRID 2 1. 0. 0. +GRID 3 -1. 0. 0. +$ELEMENTS +CBEAM 41 31 1 2 0. 1. 0. +CBEAM 42 31 1 3 0. 1. 0. +$PROPERTIES +PBEAM 31 21 10. 1. 1. 10. +$MATERIALS +MAT1 21 100. 0. +$MASSES +CONM2 11 1 5. + 5. 5. 5. +CONM2 12 2 2.5 + .00001 .00001 .00001 +CONM2 13 3 2.5 + .00001 .00001 .00001 diff --git a/examples/RafaFlyingBeam/NASTRAN/model2.bdf b/examples/RafaFlyingBeam/NASTRAN/model2.bdf new file mode 100644 index 0000000..46f6f07 --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/model2.bdf @@ -0,0 +1,23 @@ +$pyNastran: version=msc +$pyNastran: punch=True +$pyNastran: encoding=utf-8 +$pyNastran: nnodes=3 +$pyNastran: nelements=2 +$NODES +GRID 1 0. 0. 0. +GRID 2 1. 0. 0. +GRID 3 -1. 0. 0. +$ELEMENTS +CBEAM 41 31 1 2 0. 1. 0. +CBEAM 42 31 1 3 0. 1. 0. +$PROPERTIES +PBEAM 31 21 10. 1. 1. 10. +$MATERIALS +MAT1 21 100. 0. +$MASSES +CONM2 11 1 5. -.1 + 5. 5. 5. +CONM2 12 2 2.5 + .00001 .00001 .00001 +CONM2 13 3 2.5 + .00001 .00001 .00001 diff --git a/examples/RafaFlyingBeam/NASTRAN/model3.bdf b/examples/RafaFlyingBeam/NASTRAN/model3.bdf new file mode 100644 index 0000000..d81d876 --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/model3.bdf @@ -0,0 +1,23 @@ +$pyNastran: version=msc +$pyNastran: punch=True +$pyNastran: encoding=utf-8 +$pyNastran: nnodes=3 +$pyNastran: nelements=2 +$NODES +GRID 1 0. 0. 0. +GRID 2 .8660254 0. .5 +GRID 3 -.866025 0. -.5 +$ELEMENTS +CBEAM 41 31 1 2 0. 1. 0. +CBEAM 42 31 1 3 0. 1. 0. +$PROPERTIES +PBEAM 31 21 10. 1. 1. 10. +$MATERIALS +MAT1 21 100. 0. +$MASSES +CONM2 11 1 5. + 5. 5. 5. +CONM2 12 2 2.5 + .00001 .00001 .00001 +CONM2 13 3 2.5 + .00001 .00001 .00001 diff --git a/examples/RafaFlyingBeam/NASTRAN/model4.bdf b/examples/RafaFlyingBeam/NASTRAN/model4.bdf new file mode 100644 index 0000000..5ebebf0 --- /dev/null +++ b/examples/RafaFlyingBeam/NASTRAN/model4.bdf @@ -0,0 +1,23 @@ +$pyNastran: version=msc +$pyNastran: punch=True +$pyNastran: encoding=utf-8 +$pyNastran: nnodes=3 +$pyNastran: nelements=2 +$NODES +GRID 1 0. 0. 0. +GRID 2 .8660254 0. .5 +GRID 3 -.866025 0. -.5 +$ELEMENTS +CBEAM 41 31 1 2 0. 1. 0. +CBEAM 42 31 1 3 0. 1. 0. +$PROPERTIES +PBEAM 31 21 10. 1. 1. 10. +$MATERIALS +MAT1 21 100. 0. +$MASSES +CONM2 11 1 5. -.1 + 5. 5. 5. +CONM2 12 2 2.5 + .00001 .00001 .00001 +CONM2 13 3 2.5 + .00001 .00001 .00001 diff --git a/examples/RafaFlyingBeam/P1_modalsolution.py b/examples/RafaFlyingBeam/P1_modalsolution.py new file mode 100644 index 0000000..a14c4bf --- /dev/null +++ b/examples/RafaFlyingBeam/P1_modalsolution.py @@ -0,0 +1,135 @@ +# [[file:modelgen.org::*Build Nastran models][Build Nastran models:1]] +from pyNastran.bdf.bdf import BDF +import numpy as np +from dataclasses import dataclass +import pathlib +import feniax.unastran.matrixbuilder as matrixbuilder +import feniax.unastran.op2reader as op2reader +from feniax.unastran.asetbuilder import BuildAsetModel + +pathlib.Path('./FEM').mkdir(parents=True, exist_ok=True) +pathlib.Path('./NASTRAN/data_out').mkdir(parents=True, exist_ok=True) +pathlib.Path('./NASTRAN/simulations_out').mkdir(parents=True, exist_ok=True) +# Build Nastran models:1 ends here + +# [[file:modelgen.org::*Build Nastran models][Build Nastran models:2]] +@dataclass +class Config: + L: float = 1. + THETA0: float = 0. + M: float = 10. + I: float = 1. + Im: float = M * L **2 /2 + OFFSET_x: float = 0. + OFFSET_z: float = 0. + E: float = 1e2 + A: float = 5. + J: float = 10 * I + omega: float = None + + def __post_init__(self): + + self.omega = (24*self.E * self.I / (self.M * self.L ** 3) )**0.5 + +def build_bdf(config: Config): + + mesh=BDF(debug=True) + ############################ + node1 = ['GRID', 1, None, 0., 0., 0., None, None, None] + node2 = ['GRID', 2, None, config.L * np.cos(config.THETA0), 0., config.L * np.sin(config.THETA0), None, None, None] + node3 = ['GRID', 3, None, -config.L * np.cos(config.THETA0), 0., -config.L * np.sin(config.THETA0), None, None, None] + mesh.add_card(node1, 'GRID') + mesh.add_card(node2, 'GRID') + mesh.add_card(node3, 'GRID') + ############################ + # CONM2=['CONM2',Eid,RefGid,0,self.inp.mass[i][k], + # self.inp.X1[i][k],self.inp.X2[i][k],self.inp.X3[i][k],None, + # self.inp.I11[i][k],self.inp.I21[i][k], self.inp.I22[i][k], + # self.inp.I31[i][k],self.inp.I32[i][k],self.inp.I33[i][k]] + conm21 = ['CONM2', 11, 1, 0, config.M / 2, + config.OFFSET_x, 0., config.OFFSET_z, None, + config.Im, 0., config.Im, 0., 0., config.Im + ] + conm22 = ['CONM2', 12, 2, 0, config.M / 4, + 0., 0., 0. , None, + 1e-5, 0., 1e-5, 0., 0., 1e-5 + ] + conm23 = ['CONM2', 13, 3, 0, config.M / 4, + 0., 0., 0. ,None, + 1e-5, 0., 1e-5, 0., 0., 1e-5 + ] + + mesh.add_card(conm21, 'CONM2') + mesh.add_card(conm22, 'CONM2') + mesh.add_card(conm23, 'CONM2') + ############################ + # mat1 = ['MAT1',id_mat,Em,None,Nu,rho1] + mat1 = ['MAT1',21, config.E, None,None,None] + mesh.add_card(mat1, 'MAT1') + ############################ + # pbeam = ['PBEAM',id_p,id_mat,Aa,I1a,I2a,I12a,Ja] + pbeam = ['PBEAM', 31, 21, config.A, config.I, config.I, 0., config.J] + mesh.add_card(pbeam, 'PBEAM') + ############################ + # cbeam=['CBEAM',EID,PID,GA,GB,X1,X2,X3] + cbeam1= ['CBEAM', 41, 31, 1, 2, 0., 1., 0.] + cbeam2= ['CBEAM', 42, 31, 1, 3, 0., 1., 0.] + mesh.add_card(cbeam1, 'CBEAM') + mesh.add_card(cbeam2, 'CBEAM') + ############################ + return mesh +# Build Nastran models:2 ends here + +# [[file:modelgen.org::*Create nastran files for FE extraction][Create nastran files for FE extraction:1]] +config1 = Config() +mesh1 = build_bdf(config1) +mesh1.write_bdf("./NASTRAN/model1.bdf", size=8, is_double=False, close=True) +# Create nastran files for FE extraction:1 ends here + +# [[file:modelgen.org::*Create nastran files for FE extraction][Create nastran files for FE extraction:1]] +config2 = Config(OFFSET_z = -0.1) +mesh2 = build_bdf(config2) +mesh2.write_bdf("./NASTRAN/model2.bdf", size=8, is_double=False, close=True) +# Create nastran files for FE extraction:1 ends here + +# [[file:modelgen.org::*Create nastran files for FE extraction][Create nastran files for FE extraction:1]] +config3 = Config(THETA0=30*np.pi/180) +mesh3 = build_bdf(config3) +mesh3.write_bdf("./NASTRAN/model3.bdf", size=8, is_double=False, close=True) +# Create nastran files for FE extraction:1 ends here + +# [[file:modelgen.org::*Create nastran files for FE extraction][Create nastran files for FE extraction:1]] +config4 = Config(OFFSET_z = -0.1, THETA0=30*np.pi/180,) +mesh4 = build_bdf(config4) +mesh4.write_bdf("./NASTRAN/model4.bdf", size=8, is_double=False, close=True) +# Create nastran files for FE extraction:1 ends here + +# [[file:modelgen.org::*Read and save FEM and FENIAX grid][Read and save FEM and FENIAX grid:1]] +num_models = 4 +eigenvalues_list = [] +eigenvectors_list = [] +for i in range(1, num_models + 1): + op2 = op2reader.NastranReader(op2name=f"./NASTRAN/simulations_out/Model{i}_103op2.op2") + op2.readModel() + eigenvalues = op2.eigenvalues() + eigenvectors = op2.eigenvectors() + eigenvalues_list.append(eigenvalues) + eigenvectors_list.append(eigenvectors) + v = eigenvectors.reshape((18,18)).T + np.save(f"./FEM/eigenvals_m{i}.npy", eigenvalues) + np.save(f"./FEM/eigenvecs_m{i}.npy", v) + + id_list,stiffnessMatrix,massMatrix = matrixbuilder.read_pch(f"./NASTRAN/simulations_out/Model{i}_103pch.pch") + np.save(f"./FEM/Ka_m{i}.npy", stiffnessMatrix) + np.save(f"./FEM/Ma_m{i}.npy", massMatrix) +# Read and save FEM and FENIAX grid:1 ends here + +# [[file:modelgen.org::*Read and save FEM and FENIAX grid][Read and save FEM and FENIAX grid:2]] +for i in range(1, num_models + 1): + + bdf = BDF() + bdf.read_bdf(f"./NASTRAN/Model{i}_103op2.bdf", validate=False) + components = dict(rbeam=[1,2], lbeam=[3]) + model = BuildAsetModel(components, bdf) + model.write_grid(f"./FEM/structuralGrid_m{i}") +# Read and save FEM and FENIAX grid:2 ends here diff --git a/examples/RafaFlyingBeam/P2_runmodal.sh b/examples/RafaFlyingBeam/P2_runmodal.sh new file mode 100644 index 0000000..6232a11 --- /dev/null +++ b/examples/RafaFlyingBeam/P2_runmodal.sh @@ -0,0 +1,39 @@ +# [[file:modelgen.org::*Run nastran][Run nastran:1]] +source ../../feniax/unastran/run_nastran.sh +cd ./NASTRAN +run_nastran Model1_103op2.bdf +move_outputs Model1_103op2.bdf +run_nastran Model1_103pch.bdf +move_outputs Model1_103pch.bdf +cd - +# Run nastran:1 ends here + +# [[file:modelgen.org::*Run nastran][Run nastran:1]] +source ../../feniax/unastran/run_nastran.sh +cd ./NASTRAN +run_nastran Model2_103op2.bdf +move_outputs Model2_103op2.bdf +run_nastran Model2_103pch.bdf +move_outputs Model2_103pch.bdf +cd - +# Run nastran:1 ends here + +# [[file:modelgen.org::*Run nastran][Run nastran:1]] +source ../../feniax/unastran/run_nastran.sh +cd ./NASTRAN +run_nastran Model3_103op2.bdf +move_outputs Model3_103op2.bdf +run_nastran Model3_103pch.bdf +move_outputs Model3_103pch.bdf +cd - +# Run nastran:1 ends here + +# [[file:modelgen.org::*Run nastran][Run nastran:1]] +source ../../feniax/unastran/run_nastran.sh +cd ./NASTRAN +run_nastran Model4_103op2.bdf +move_outputs Model4_103op2.bdf +run_nastran Model4_103pch.bdf +move_outputs Model4_103pch.bdf +cd - +# Run nastran:1 ends here diff --git a/examples/RafaFlyingBeam/P3_simulation.py b/examples/RafaFlyingBeam/P3_simulation.py new file mode 100644 index 0000000..c3b6477 --- /dev/null +++ b/examples/RafaFlyingBeam/P3_simulation.py @@ -0,0 +1,53 @@ +# [[file:modelgen.org::*FENIAX][FENIAX:1]] +import feniax.preprocessor.configuration as configuration +from feniax.preprocessor.inputs import Inputs +import feniax.feniax_main +import jax.numpy as jnp +import pathlib +# FENIAX:1 ends here + +# [[file:modelgen.org::*FENIAX][FENIAX:2]] +v_x = 0. +v_y = 0. +v_z = 0. +omega_x = 0. +omega_y = 1. +omega_z = 0. +gravity_forces = False +label = 'm1' +# FENIAX:2 ends here + +# [[file:modelgen.org::*FENIAX][FENIAX:3]] +inp = Inputs() +inp.engine = "intrinsicmodal" +inp.fem.connectivity = {'rbeam': None, 'lbeam': None} +inp.fem.Ka_name = f"./FEM/Ka_{label}.npy" +inp.fem.Ma_name = f"./FEM/Ma_{label}.npy" +inp.fem.eig_names = [f"./FEM/eigenvals_{label}.npy", + f"./FEM/eigenvecs_{label}.npy"] +inp.fem.grid = f"./FEM/structuralGrid_{label}" +inp.fem.num_modes = 18 +inp.fem.eig_type = "inputs" +inp.driver.typeof = "intrinsic" +inp.driver.sol_path= pathlib.Path( + f"./results_{label}") +inp.simulation.typeof = "single" +inp.system.name = "s1" +inp.system.solution = "dynamic" +inp.system.bc1 = 'free' +inp.system.xloads.gravity_forces = gravity_forces +inp.system.t1 = 2. +inp.system.tn = 2001 +inp.system.solver_library = "runge_kutta" #"diffrax" # +inp.system.solver_function = "ode" +inp.system.solver_settings = dict(solver_name="rk4") +inp.system.init_states = dict(q1=["nodal_prescribed", + ([[v_x, v_y, v_z, omega_x, omega_y, omega_z], + [v_x, v_y, v_z, omega_x, omega_y, omega_z], + [v_x, v_y, v_z, omega_x, omega_y, omega_z]] + ,) + ] + ) +config = configuration.Config(inp) +sol = feniax.feniax_main.main(input_obj=config) +# FENIAX:3 ends here diff --git a/examples/RafaFlyingBeam/Streamlit/Home.py b/examples/RafaFlyingBeam/Streamlit/Home.py new file mode 100644 index 0000000..2d27894 --- /dev/null +++ b/examples/RafaFlyingBeam/Streamlit/Home.py @@ -0,0 +1,38 @@ +#Import the required Libraries +import streamlit as st +from feniax.preprocessor import solution +import feniax.preprocessor.configuration as configuration +import importlib +import feniax.plotools.streamlit.intrinsic as sti +importlib.reload(sti) +import feniax + +st.set_page_config( + page_title="Home page", + page_icon="🖥️", + layout="wide" + ) + +sti.home() +st.divider() +st.title("") +# st.markdown("This is a simplified box wing made of composite shells, lumped masses along the load paths and interpolation elements (RBE3s) connecting those to the main structure.") +# wingsp_path = str(feniax.PATH / "../docs/images/wingSP5b.png") +# st.image(wingsp_path, caption="Wing box FE model") + + +st.text("This is what we can do!!") + +st.markdown(""" +### Select a folder with results for postprocessing +""") +solfolder = "../results_m1" +selected_folder = sti.file_selector('../') +if selected_folder is not None: + solfolder = selected_folder +st.write('Solution Folder `%s`' % solfolder) +if solfolder is not None: + sol = solution.IntrinsicReader(solfolder) + config = configuration.Config.from_file(f"{solfolder}/config.yaml") + st.session_state['sol'] = sol + st.session_state['config'] = config diff --git a/examples/RafaFlyingBeam/Streamlit/pages/1_Geometry.py b/examples/RafaFlyingBeam/Streamlit/pages/1_Geometry.py new file mode 100644 index 0000000..1cddc8a --- /dev/null +++ b/examples/RafaFlyingBeam/Streamlit/pages/1_Geometry.py @@ -0,0 +1,25 @@ +import feniax.plotools.streamlit.intrinsic as sti +import streamlit as st +import feniax.plotools.streamlit.theory as stt + +import importlib +importlib.reload(sti) + +st.set_page_config( + page_title="Initial model geometry", + page_icon="🛫", + layout="wide" +) + +stt.fe_reduction() + +st.divider() +sti.df_geometry(st.session_state.config.fem) +st.divider() +sti.sys_3Dconfiguration0(st.session_state.config) +st.divider() +left, = st.columns(1) +if left.button("Click to see FE matrices", use_container_width=True): + + sti.fe_matrices(st.session_state.config.fem) + diff --git a/examples/RafaFlyingBeam/Streamlit/pages/2_Modes.py b/examples/RafaFlyingBeam/Streamlit/pages/2_Modes.py new file mode 100644 index 0000000..7bb5e7f --- /dev/null +++ b/examples/RafaFlyingBeam/Streamlit/pages/2_Modes.py @@ -0,0 +1,20 @@ +import feniax.plotools.streamlit.intrinsic as sti +import feniax.plotools.streamlit.theory as stt + +import streamlit as st +import importlib +importlib.reload(sti) + +st.set_page_config( + page_title="Intrinsic modal shapes", + page_icon="🛫", + layout="wide" +) + +stt.intrinsic_modes() +st.divider() +st.header('Intrinsic modal data') +st.link_button("Code","https://github.com/ACea15/FENIAX/blob/0958ca92b55073d799668136ae4d5132687f8969/feniax.intrinsic/modes.py#L192") + +sti.df_modes(st.session_state.sol, + st.session_state.config) diff --git a/examples/RafaFlyingBeam/Streamlit/pages/3_ModalCouplings.py b/examples/RafaFlyingBeam/Streamlit/pages/3_ModalCouplings.py new file mode 100644 index 0000000..b334773 --- /dev/null +++ b/examples/RafaFlyingBeam/Streamlit/pages/3_ModalCouplings.py @@ -0,0 +1,21 @@ +import feniax.plotools.streamlit.intrinsic as sti +import streamlit as st +import importlib +importlib.reload(sti) +import feniax.plotools.streamlit.theory as stt + +st.set_page_config( + page_title="Intrinsic modal couplings", + page_icon="🚁", + layout="wide" +) + +stt.intrinsic_couplings() +st.divider() +st.header('Modal couplings') + +st.link_button("Code","https://github.com/ACea15/FENIAX/blob/a54b758c10b53e203268a810d6bf813160b34320/fem4inas/intrinsic/couplings.py#L9") + + +st.divider() +sti.df_couplings(st.session_state.sol) diff --git a/examples/RafaFlyingBeam/Streamlit/pages/4_Systems.py b/examples/RafaFlyingBeam/Streamlit/pages/4_Systems.py new file mode 100644 index 0000000..812ca1d --- /dev/null +++ b/examples/RafaFlyingBeam/Streamlit/pages/4_Systems.py @@ -0,0 +1,15 @@ +import feniax.plotools.streamlit.intrinsic as sti +import streamlit as st +# import importlib +# importlib.reload(sti) +import feniax.plotools.streamlit.theory as stt + +st.set_page_config( + page_title="Systems", + page_icon="🚁", + layout="wide" +) + +stt.intrinsic_systems() + +sti.systems(st.session_state.sol, st.session_state.config) diff --git a/examples/RafaFlyingBeam/Streamlit/pages/5_Comparison.py b/examples/RafaFlyingBeam/Streamlit/pages/5_Comparison.py new file mode 100644 index 0000000..8a5d08c --- /dev/null +++ b/examples/RafaFlyingBeam/Streamlit/pages/5_Comparison.py @@ -0,0 +1,34 @@ +import feniax.plotools.streamlit.intrinsic as sti +from feniax.preprocessor import solution +import feniax.preprocessor.configuration as configuration + +import streamlit as st +import importlib +importlib.reload(sti) + +st.set_page_config( + page_title="Comparison", + page_icon="🚁", + layout="wide" +) + + +Csol = dict() +Cconfig = dict() +st.markdown(""" +### Select folders with results for postprocessing +""") +selected_folders = sti.multifolder_selector('../') +if len(selected_folders) > 0: + #breakpoint() + for si in selected_folders: + #st.write('Solution Folder `%s`' % solfolder) + sol = solution.IntrinsicReader(si) + config = configuration.Config.from_file(f"{si}/config.yaml") + Csol[si.name] = sol + Cconfig[si.name] = config + st.session_state['Csol'] = Csol + st.session_state['Cconfig'] = Cconfig + sti.systems_comparison(st.session_state['Csol'], + st.session_state['Cconfig']) +#sti.systems(st.session_state.sol, st.session_state.config) diff --git a/examples/RafaFlyingBeam/Streamlit/pages/6_Background.py b/examples/RafaFlyingBeam/Streamlit/pages/6_Background.py new file mode 100644 index 0000000..8773218 --- /dev/null +++ b/examples/RafaFlyingBeam/Streamlit/pages/6_Background.py @@ -0,0 +1,32 @@ +#Import the required Libraries +import streamlit as st +from feniax.preprocessor import solution +import feniax.preprocessor.configuration as configuration +import importlib +import feniax.plotools.streamlit.intrinsic as sti +importlib.reload(sti) +import fem4inas + +st.set_page_config( + page_title="Background page", + page_icon="🖥️", + layout="wide" + ) + +st.title("Background theory") + +st.header("Solution process") +path = feniax.PATH / "../docs/images" +st.image(str(path / "aircraft_process2.png"), + caption="Steps for the structural solution") +st.divider() + +st.header("Theoretical assumptions of the solution") +st.image(str(path / "reality2NMROM2.png"), + caption="Underlying assumptions in the solution") +st.divider() + +st.header("Implementation in JAX") +st.image(str(path / "jaxlogo.png"), + caption="JAX capabilities") +st.divider() diff --git a/examples/RafaFlyingBeam/img/3mass_system.png b/examples/RafaFlyingBeam/img/3mass_system.png new file mode 100644 index 0000000..463e8e9 Binary files /dev/null and b/examples/RafaFlyingBeam/img/3mass_system.png differ diff --git a/examples/RafaFlyingBeam/modelgen.org b/examples/RafaFlyingBeam/modelgen.org new file mode 100644 index 0000000..3325d98 --- /dev/null +++ b/examples/RafaFlyingBeam/modelgen.org @@ -0,0 +1,377 @@ +#+TITLE: Example 6.1 of Dynamics of Flexible Aircraft + +* House keeping +#+begin_src elisp :results none + (setq org-confirm-babel-evaluate nil) + (require 'org-tempo) + (pyvenv-workon "feniax") +#+end_src + +* Problem description + +Problem: free dynamics of a 3-mass flexible system + +[[./img/3mass_system.png]] + + +Parameters in the problem: +- L: Arm +- THETA0: Initial inclination +- M: parameter mass +- I: Moment of Inertia (Izz & Iyy) +- J: Moment of inertia Ixx +- Im: Mass moment of inertia = M * L **2 /2 +- OFFSET_x: Offset in x of the central mass with respect to the grid node +- OFFSET_z: Offset in z of the central mass with respect to the grid node +- E: young modulus +- A: cross-sectional area +- omega: theoretical natural frquency + + We will build various FE models in Nastran first, then simulate the free dynamics in FENIAX. + +* Build Nastran models +:PROPERTIES: +:header-args: :tangle ./P1_modalsolution.py :session *pyshell* :comments yes :results none +:END: + +#+begin_src python + from pyNastran.bdf.bdf import BDF + import numpy as np + from dataclasses import dataclass + import pathlib + import feniax.unastran.matrixbuilder as matrixbuilder + import feniax.unastran.op2reader as op2reader + from feniax.unastran.asetbuilder import BuildAsetModel + + pathlib.Path('./FEM').mkdir(parents=True, exist_ok=True) + pathlib.Path('./NASTRAN/data_out').mkdir(parents=True, exist_ok=True) + pathlib.Path('./NASTRAN/simulations_out').mkdir(parents=True, exist_ok=True) + +#+end_src + + +#+begin_src python + + @dataclass + class Config: + L: float = 1. + THETA0: float = 0. + M: float = 10. + I: float = 1. + Im: float = M * L **2 /2 + OFFSET_x: float = 0. + OFFSET_z: float = 0. + E: float = 1e2 + A: float = 5. + J: float = 10 * I + omega: float = None + + def __post_init__(self): + + self.omega = (24*self.E * self.I / (self.M * self.L ** 3) )**0.5 + + def build_bdf(config: Config): + + mesh=BDF(debug=True) + ############################ + node1 = ['GRID', 1, None, 0., 0., 0., None, None, None] + node2 = ['GRID', 2, None, config.L * np.cos(config.THETA0), 0., config.L * np.sin(config.THETA0), None, None, None] + node3 = ['GRID', 3, None, -config.L * np.cos(config.THETA0), 0., -config.L * np.sin(config.THETA0), None, None, None] + mesh.add_card(node1, 'GRID') + mesh.add_card(node2, 'GRID') + mesh.add_card(node3, 'GRID') + ############################ + # CONM2=['CONM2',Eid,RefGid,0,self.inp.mass[i][k], + # self.inp.X1[i][k],self.inp.X2[i][k],self.inp.X3[i][k],None, + # self.inp.I11[i][k],self.inp.I21[i][k], self.inp.I22[i][k], + # self.inp.I31[i][k],self.inp.I32[i][k],self.inp.I33[i][k]] + conm21 = ['CONM2', 11, 1, 0, config.M / 2, + config.OFFSET_x, 0., config.OFFSET_z, None, + config.Im, 0., config.Im, 0., 0., config.Im + ] + conm22 = ['CONM2', 12, 2, 0, config.M / 4, + 0., 0., 0. , None, + 1e-5, 0., 1e-5, 0., 0., 1e-5 + ] + conm23 = ['CONM2', 13, 3, 0, config.M / 4, + 0., 0., 0. ,None, + 1e-5, 0., 1e-5, 0., 0., 1e-5 + ] + + mesh.add_card(conm21, 'CONM2') + mesh.add_card(conm22, 'CONM2') + mesh.add_card(conm23, 'CONM2') + ############################ + # mat1 = ['MAT1',id_mat,Em,None,Nu,rho1] + mat1 = ['MAT1',21, config.E, None,None,None] + mesh.add_card(mat1, 'MAT1') + ############################ + # pbeam = ['PBEAM',id_p,id_mat,Aa,I1a,I2a,I12a,Ja] + pbeam = ['PBEAM', 31, 21, config.A, config.I, config.I, 0., config.J] + mesh.add_card(pbeam, 'PBEAM') + ############################ + # cbeam=['CBEAM',EID,PID,GA,GB,X1,X2,X3] + cbeam1= ['CBEAM', 41, 31, 1, 2, 0., 1., 0.] + cbeam2= ['CBEAM', 42, 31, 1, 3, 0., 1., 0.] + mesh.add_card(cbeam1, 'CBEAM') + mesh.add_card(cbeam2, 'CBEAM') + ############################ + return mesh + +#+end_src + +#+NAME: bdf103bulk +#+begin_src org :tangle no + SOL 103 + CEND + TITLE=Rafa's flying beam model # + ECHO=NONE + DISPLACEMENT=ALL + METHOD = 900 + BEGIN BULK + EIGRL,900,,,18 +#+end_src + +** Model 1 + +Horizontal bar, no offset +*** Create nastran files for FE extraction +#+begin_src python + config1 = Config() + mesh1 = build_bdf(config1) + mesh1.write_bdf("./NASTRAN/model1.bdf", size=8, is_double=False, close=True) +#+end_src + +- For eigenvectors: +#+begin_src org :noweb yes :tangle ./NASTRAN/Model1_103op2.bdf :comments no + <> + PARAM,POST,-1 + INCLUDE 'model1.bdf' + ENDDATA +#+end_src + +- pch for FE matrices +#+begin_src org :noweb yes :tangle ./NASTRAN/Model1_103pch.bdf :comments no + <> + PARAM,EXTOUT,DMIGPCH + INCLUDE 'model1.bdf' + ENDDATA +#+end_src + +*** Run nastran +#+begin_src bash :session shell1 :tangle P2_runmodal.sh + source ../../feniax/unastran/run_nastran.sh + cd ./NASTRAN + run_nastran Model1_103op2.bdf + move_outputs Model1_103op2.bdf + run_nastran Model1_103pch.bdf + move_outputs Model1_103pch.bdf + cd - +#+end_src + +** Model 2 + +Horizontal bar, 0.1 offset +*** Create nastran files for FE extraction +#+begin_src python + config2 = Config(OFFSET_z = -0.1) + mesh2 = build_bdf(config2) + mesh2.write_bdf("./NASTRAN/model2.bdf", size=8, is_double=False, close=True) +#+end_src + +- For eigenvectors: +#+begin_src org :noweb yes :tangle ./NASTRAN/Model2_103op2.bdf :comments no + <> + PARAM,POST,-1 + INCLUDE 'model2.bdf' + ENDDATA +#+end_src + +- pch for FE matrices +#+begin_src org :noweb yes :tangle ./NASTRAN/Model2_103pch.bdf :comments no + <> + PARAM,EXTOUT,DMIGPCH + INCLUDE 'model2.bdf' + ENDDATA +#+end_src + +*** Run nastran +#+begin_src bash :session shell1 :tangle P2_runmodal.sh + source ../../feniax/unastran/run_nastran.sh + cd ./NASTRAN + run_nastran Model2_103op2.bdf + move_outputs Model2_103op2.bdf + run_nastran Model2_103pch.bdf + move_outputs Model2_103pch.bdf + cd - +#+end_src + +** Model 3 + +Inclined 30 degrees bar, no offset +*** Create nastran files for FE extraction +#+begin_src python + config3 = Config(THETA0=30*np.pi/180) + mesh3 = build_bdf(config3) + mesh3.write_bdf("./NASTRAN/model3.bdf", size=8, is_double=False, close=True) +#+end_src + +- For eigenvectors: +#+begin_src org :noweb yes :tangle ./NASTRAN/Model3_103op2.bdf :comments no + <> + PARAM,POST,-1 + INCLUDE 'model3.bdf' + ENDDATA +#+end_src + +- pch for FE matrices +#+begin_src org :noweb yes :tangle ./NASTRAN/Model3_103pch.bdf :comments no + <> + PARAM,EXTOUT,DMIGPCH + INCLUDE 'model3.bdf' + ENDDATA +#+end_src + +*** Run nastran +#+begin_src bash :session shell1 :tangle P2_runmodal.sh + source ../../feniax/unastran/run_nastran.sh + cd ./NASTRAN + run_nastran Model3_103op2.bdf + move_outputs Model3_103op2.bdf + run_nastran Model3_103pch.bdf + move_outputs Model3_103pch.bdf + cd - +#+end_src + +** Model 4 +Inclined 30 degrees bar, 0.1 offset +*** Create nastran files for FE extraction +#+begin_src python + config4 = Config(OFFSET_z = -0.1, THETA0=30*np.pi/180,) + mesh4 = build_bdf(config4) + mesh4.write_bdf("./NASTRAN/model4.bdf", size=8, is_double=False, close=True) +#+end_src + +- For eigenvectors: +#+begin_src org :noweb yes :tangle ./NASTRAN/Model4_103op2.bdf :comments no + <> + PARAM,POST,-1 + INCLUDE 'model4.bdf' + ENDDATA +#+end_src + +- pch for FE matrices +#+begin_src org :noweb yes :tangle ./NASTRAN/Model4_103pch.bdf :comments no + <> + PARAM,EXTOUT,DMIGPCH + INCLUDE 'model4.bdf' + ENDDATA +#+end_src + +*** Run nastran +#+begin_src bash :session shell1 :tangle P2_runmodal.sh + source ../../feniax/unastran/run_nastran.sh + cd ./NASTRAN + run_nastran Model4_103op2.bdf + move_outputs Model4_103op2.bdf + run_nastran Model4_103pch.bdf + move_outputs Model4_103pch.bdf + cd - +#+end_src + +** Read and save FEM and FENIAX grid + +- Save Ka, Ma, eigenvalues and eigenvectors +#+begin_src python + num_models = 4 + eigenvalues_list = [] + eigenvectors_list = [] + for i in range(1, num_models + 1): + op2 = op2reader.NastranReader(op2name=f"./NASTRAN/simulations_out/Model{i}_103op2.op2") + op2.readModel() + eigenvalues = op2.eigenvalues() + eigenvectors = op2.eigenvectors() + eigenvalues_list.append(eigenvalues) + eigenvectors_list.append(eigenvectors) + v = eigenvectors.reshape((18,18)).T + np.save(f"./FEM/eigenvals_m{i}.npy", eigenvalues) + np.save(f"./FEM/eigenvecs_m{i}.npy", v) + + id_list,stiffnessMatrix,massMatrix = matrixbuilder.read_pch(f"./NASTRAN/simulations_out/Model{i}_103pch.pch") + np.save(f"./FEM/Ka_m{i}.npy", stiffnessMatrix) + np.save(f"./FEM/Ma_m{i}.npy", massMatrix) +#+end_src + +- save Grid file +#+begin_src python + for i in range(1, num_models + 1): + + bdf = BDF() + bdf.read_bdf(f"./NASTRAN/Model{i}_103op2.bdf", validate=False) + components = dict(rbeam=[1,2], lbeam=[3]) + model = BuildAsetModel(components, bdf) + model.write_grid(f"./FEM/structuralGrid_m{i}") + +#+end_src + +* FENIAX +:PROPERTIES: +:header-args: :tangle ./P3_simulation.py :session *pyfeniax* :comments yes :results none +:END: + +Load simulation modules +#+begin_src python + import feniax.preprocessor.configuration as configuration + from feniax.preprocessor.inputs import Inputs + import feniax.feniax_main + import jax.numpy as jnp + import pathlib +#+end_src + +Set model to be run (mi), initial conditions and whether to include gravity forces: +#+begin_src python + v_x = 1. + v_y = 0. + v_z = 0. + omega_x = 0. + omega_y = 1. + omega_z = 0. + gravity_forces = False + label = 'm1' +#+end_src + +#+begin_src python + inp = Inputs() + inp.engine = "intrinsicmodal" + inp.fem.connectivity = {'rbeam': None, 'lbeam': None} + inp.fem.Ka_name = f"./FEM/Ka_{label}.npy" + inp.fem.Ma_name = f"./FEM/Ma_{label}.npy" + inp.fem.eig_names = [f"./FEM/eigenvals_{label}.npy", + f"./FEM/eigenvecs_{label}.npy"] + inp.fem.grid = f"./FEM/structuralGrid_{label}" + inp.fem.num_modes = 18 + inp.fem.eig_type = "inputs" + inp.driver.typeof = "intrinsic" + inp.driver.sol_path= pathlib.Path( + f"./results_{label}") + inp.simulation.typeof = "single" + inp.system.name = "s1" + inp.system.solution = "dynamic" + inp.system.bc1 = 'free' + inp.system.xloads.gravity_forces = gravity_forces + inp.system.t1 = 1. + inp.system.tn = 5001 + inp.system.solver_library = "runge_kutta" #"diffrax" # + inp.system.solver_function = "ode" + inp.system.solver_settings = dict(solver_name="rk4") + inp.system.init_states = dict(q1=["nodal_prescribed", + ([[v_x, v_y, v_z, omega_x, omega_y, omega_z], + [v_x, v_y, v_z, omega_x, omega_y, omega_z], + [v_x, v_y, v_z, omega_x, omega_y, omega_z]] + ,) + ] + ) + config = configuration.Config(inp) + sol = feniax.feniax_main.main(input_obj=config) + +#+end_src