Allowing user to Initialize pyTACS from pynastran BDF object (#245)

* Adding option to initialize pyTACS from pyNASTRAN BDF object

* Adding example to shown new intialization feature

* Adding pynastran initialization test

* Adding docstrings to example

* Minor error reformat

* Another minor error reformat

* Tightening real test tol

examples/annulus/analysis.py

@@ -0,0 +1,99 @@
+"""
+This example shows how to initialize pyTACS directly from a pyNASTRAN BDF object.
+In this example we create an annulus out of shell elements using pyNASTRAN's BDF class.
+We then use this class object to initialize pyTACS rather than reading in a file.
+We analyze two load cases: A static case where we apply a load at the spider-webbed RBE2 centered inside the annulus
+and modal case where we evaluate the mode shapes of the structure.
+"""
+
+import numpy as np
+from pyNastran.bdf.bdf import BDF
+
+from tacs import pyTACS, functions
+
+
+# Geometric properties
+Ri = 0.5
+Ro = 1.0
+z = 0.0
+t = 0.05
+# Material properties
+E = 70e9
+nu = 0.3
+ys = 270e6
+rho = 2700.0
+
+# number of elements in r/theta
+nr = 10
+ntheta = 20
+# Center node ID
+center_nid = 1000
+
+# Initialize our BDF object
+bdfInfo = BDF()
+# Create a cylindrical coord system
+bdfInfo.add_cord2c(1, [0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 1.0])
+# Discretize in r and theta
+r_array = np.linspace(Ri, Ro, nr)
+theta_array = np.linspace(0.0, 360.0, ntheta, endpoint=False)
+
+
+# Loop through and add each node
+curID = 1
+nodeIDs = {}
+for r in r_array:
+    for theta in theta_array:
+        bdfInfo.add_grid(curID, [r, theta, z], 1)
+        nodeIDs[r, theta] = curID
+        curID += 1
+
+# Loop through and add each element
+curID = 1
+for i in range(nr - 1):
+    for j in range(ntheta):
+        r0 = r_array[i]
+        t0 = theta_array[j - 1]
+        r1 = r_array[i + 1]
+        t1 = theta_array[j]
+        conn = [nodeIDs[r0, t0], nodeIDs[r0, t1], nodeIDs[r1, t1], nodeIDs[r1, t0]]
+        bdfInfo.add_cquad4(curID, 1, conn)
+        curID += 1
+
+# Simply support the outer edge
+for theta in theta_array:
+    bdfInfo.add_spc(0, nodeIDs[Ro, theta], "123", 0.0)
+
+# Add node at center for applying load
+bdfInfo.add_grid(center_nid, [0.0, 0.0, 0.0], 0)
+
+# Connect center node with RBE to inner edge
+rbe_dep_nodes = []
+for theta in theta_array:
+    rbe_dep_nodes.append(nodeIDs[Ri, theta])
+bdfInfo.add_rbe2(curID, center_nid, "123456", rbe_dep_nodes)
+
+# Define material card for Aluminum
+bdfInfo.add_mat1(1, E, None, nu, rho, St=ys)
+# Define shell thickness
+bdfInfo.add_pshell(1, 1, t)
+
+# Initialize pyTACS from our BDF object
+FEAAssembler = pyTACS(bdfInfo)
+FEAAssembler.initialize()
+
+# Apply a static load case with a x moment at the center
+prob = FEAAssembler.createStaticProblem("center_moment")
+prob.addLoadToNodes(center_nid, [0.0, 0.0, 0.0, 1e7, 0.0, 0.0], nastranOrdering=True)
+prob.addFunction(
+    "Izz", functions.MomentOfInertia, direction1=[0, 0, 1], direction2=[0, 0, 1]
+)
+prob.solve()
+f = {}
+prob.evalFunctions(f)
+print(f)
+prob.writeSolution()
+
+# Perform modal analysis
+modal = FEAAssembler.createModalProblem("modes", 1.0, 10)
+modal.solve()
+modal.writeSolution()

tacs/pymeshloader.py

@@ -10,10 +10,11 @@
 # =============================================================================
 # Imports
 # =============================================================================
+from copy import deepcopy
 import itertools as it
 
 import numpy as np
-from pyNastran.bdf.bdf import read_bdf
+import pyNastran.bdf.bdf as pn
 
 import tacs.TACS
 import tacs.constitutive
@@ -29,10 +30,10 @@ def __init__(self, comm, printDebug=False):
         self.printDebug = printDebug
         self.bdfInfo = None
 
-    def scanBdfFile(self, fileName):
+    def scanBdfFile(self, bdf):
         """
         Scan nastran bdf file using pyNastran's bdf parser.
-        We also set up arrays that will be require later to build tacs.
+        We also set up arrays that will be required later to build tacs.
         """
 
         # Only print debug info on root, if requested
@@ -41,10 +42,20 @@ def scanBdfFile(self, fileName):
         else:
             debugPrint = False
 
-        # Read in bdf file as pynastran object
-        # By default we avoid cross-referencing unless we actually need it,
-        # since its expensive for large models
-        self.bdfInfo = read_bdf(fileName, validate=False, xref=False, debug=debugPrint)
+        # Check if a file name was provided
+        # Create a BDF object from the file
+        if isinstance(bdf, str):
+            # Read in bdf file as pynastran object
+            # By default we avoid cross-referencing unless we actually need it,
+            # since its expensive for large models
+            self.bdfInfo = pn.read_bdf(bdf, validate=False, xref=False, debug=debugPrint)
+        # Create a copy of the BDF object
+        elif isinstance(bdf, pn.BDF):
+            self.bdfInfo = deepcopy(bdf)
+        else:
+            raise self._TACSError("BDF input must be provided as a file name 'str' or pyNastran 'BDF' object. "
+                                  f"Provided input was of type '{type(bdf).__name__}'.")
+
         # Set flag letting us know model is not xrefed yet
         self.bdfInfo.is_xrefed = False
 

tacs/pytacs.py

@@ -153,13 +153,13 @@ class pyTACS(BaseUI):
         ],
     }
 
-    def __init__(self, fileName, comm=None, dvNum=0, scaleList=None, options=None):
+    def __init__(self, bdf, comm=None, dvNum=0, scaleList=None, options=None):
         """
 
         Parameters
         ----------
-        fileName : str
-            The filename of the BDF file to load.
+        bdf : str or pyNastran.bdf.bdf.BDF
+            The BDF file or a pyNastran BDF object to load.
 
         comm : mpi4py.MPI.Intracomm
             The comm object on which to create the pyTACS object.
@@ -189,10 +189,10 @@ def __init__(self, fileName, comm=None, dvNum=0, scaleList=None, options=None):
         # Create and load mesh loader object.
         debugFlag = self.getOption("printDebug")
         self.meshLoader = pyMeshLoader(self.comm, debugFlag)
-        self.meshLoader.scanBdfFile(fileName)
-        self.bdfName = fileName
+        self.meshLoader.scanBdfFile(bdf)
         # Save pynastran bdf object
         self.bdfInfo = self.meshLoader.getBDFInfo()
+        self.bdfName = self.bdfInfo.bdf_filename
 
         meshLoadTime = time.time()
 

tests/integration_tests/test_shell_annulus_from_pynastran.py

@@ -0,0 +1,139 @@
+import numpy as np
+from pyNastran.bdf.bdf import BDF
+
+from pytacs_analysis_base_test import PyTACSTestCase
+from tacs import pytacs, functions
+
+"""
+Create annulus out of shell elements using pyNASTRAN's BDF class.
+We then use this class object to initialize pyTACS rather than reading in a file.
+Apply a load at the spider-webbed RBE2 centered inside the annulus and test KSFailure, StructuralMass, 
+Moment of Inertia, and Compliance functions and sensitivities.
+"""
+
+# Geometric properties
+Ri = 0.5
+Ro = 1.0
+z = 0.0
+t = 0.05
+# Material properties
+E = 70e9
+nu = 0.3
+ys = 270e6
+rho = 2700.0
+
+# number of elements in r/theta
+nr = 5
+ntheta = 10
+# Center node ID
+center_nid = 1000
+
+# KS function weight
+ksweight = 10.0
+
+
+class ProblemTest(PyTACSTestCase.PyTACSTest):
+    N_PROCS = 2  # this is how many MPI processes to use for this TestCase.
+
+    FUNC_REFS = {
+        "center_moment_Izz": 174.1393226319456,
+        "center_moment_compliance": 2047161.6725086374,
+        "center_moment_ks_vmfailure": 15.179071361242388,
+        "center_moment_mass": 297.56628397306457,
+    }
+
+    def setup_tacs_problems(self, comm):
+        """
+        Setup pytacs object for problems we will be testing.
+        """
+
+        # Overwrite default check values
+        if self.dtype == complex:
+            self.rtol = 5e-6
+            self.atol = 1e-8
+            self.dh = 1e-10
+        else:
+            self.rtol = 1e-3
+            self.atol = 1e-3
+            self.dh = 1e-6
+
+        # Instantiate FEA Assembler
+        bdf_info = self.createBDFObject()
+
+        fea_assembler = pytacs.pyTACS(bdf_info, comm)
+
+        # Set up constitutive objects and elements
+        fea_assembler.initialize()
+
+        # Apply a static load case with a x moment at the center
+        problem = fea_assembler.createStaticProblem("center_moment")
+        problem.addLoadToNodes(
+            center_nid, [0.0, 0.0, 0.0, 1e7, 0.0, 0.0], nastranOrdering=True
+        )
+
+        problem.addFunction("mass", functions.StructuralMass)
+        problem.addFunction(
+            "Izz", functions.MomentOfInertia, direction1=[0, 0, 1], direction2=[0, 0, 1]
+        )
+        problem.addFunction("ks_vmfailure", functions.KSFailure, ksWeight=ksweight)
+        problem.addFunction("compliance", functions.Compliance)
+        # Set convergence tol to be tight
+        problem.setOption("L2Convergence", 1e-20)
+        problem.setOption("L2ConvergenceRel", 1e-20)
+
+        return [problem], fea_assembler
+
+    def createBDFObject(self):
+        # Initialize our BDF object
+        bdfInfo = BDF()
+        # Create a cylindrical coord system
+        bdfInfo.add_cord2c(1, [0.0, 0.0, 0.0], [0.0, 0.0, 1.0], [1.0, 0.0, 1.0])
+        # Discretize in r and theta
+        r_array = np.linspace(Ri, Ro, nr)
+        theta_array = np.linspace(0.0, 360.0, ntheta, endpoint=False)
+
+        # Loop through and add each node
+        curID = 1
+        nodeIDs = {}
+        for r in r_array:
+            for theta in theta_array:
+                bdfInfo.add_grid(curID, [r, theta, z], 1)
+                nodeIDs[r, theta] = curID
+                curID += 1
+
+        # Loop through and add each element
+        curID = 1
+        for i in range(nr - 1):
+            for j in range(ntheta):
+                r0 = r_array[i]
+                t0 = theta_array[j - 1]
+                r1 = r_array[i + 1]
+                t1 = theta_array[j]
+                conn = [
+                    nodeIDs[r0, t0],
+                    nodeIDs[r0, t1],
+                    nodeIDs[r1, t1],
+                    nodeIDs[r1, t0],
+                ]
+                bdfInfo.add_cquad4(curID, 1, conn)
+                curID += 1
+
+        # Simply support the outer edge
+        for theta in theta_array:
+            bdfInfo.add_spc(0, nodeIDs[Ro, theta], "123", 0.0)
+
+        # Add node at center for applying load
+        bdfInfo.add_grid(center_nid, [0.0, 0.0, 0.0], 0)
+
+        # Connect center node with RBE to inner edge
+        rbe_dep_nodes = []
+        for theta in theta_array:
+            rbe_dep_nodes.append(nodeIDs[Ri, theta])
+        bdfInfo.add_rbe2(curID, center_nid, "123456", rbe_dep_nodes)
+
+        # Define material card for Aluminum
+        bdfInfo.add_mat1(1, E, None, nu, rho, St=ys)
+        # Define shell thickness
+        bdfInfo.add_pshell(1, 1, t)
+
+        return bdfInfo