Adding BDF writer to pyTACS (smdogroup#198)

* Adding BDF writing methods for MaterialProperties class

* Adding BDF writing methods for IsoShellConstitutive class

* Adding BDF writing methods for SolidConstitutive class

* Adding BDF writing methods for OrthotropicPly class

* Adding BDF writing methods for CompositeShellConstitutive

* Adding BDF writing methods for BasicBeamConstitutive

* Adding BDF writing methods for DOFSpringConstitutive

* Adding getConstitutive

* Adding writeBDF method to pyTACS

* Adding minor comments

* Adding concentrated masses to bdf writer

* Adding design variable update to bdf writer

* Removing writebdf call in crm example

* Adding more comments

* Adding some docstrings to cython

* Adding PROD to cython

* Minor reodering of `writeBDF`

* Fixing comment handling in bdfwriter

* Adding load writer to StaticProblem class

* Making bdf writing methods complex-safe

* Adding boundary conditions to bdf header

* Adding new test for writeBDF

* Modifying beam procedure in writeBDF

* Adding new bdf writer test

* Adding new bdf writer to mphys

* Adding more docstrings to cython

* Minor docstring fix

* Minor docstring fix

* Updating mass coments in `writeBDF`

* Adding composite test to integration tests

* Adding composite writebdf test to integration tests

* Minor edits to `writeBDF` procedure

* Tweaking composite test

* Moving mphys write_bdf function to coupling group

examples/beam/beam_opt.py

@@ -130,6 +130,10 @@ def setup(self):
 # Run optimization
 prob.run_driver()
 
+# Write optimized structure to BDF
+bdf_out = os.path.join(os.path.dirname(__file__), "beam_sol.bdf")
+prob.model.tip_shear.coupling.write_bdf(bdf_out)
+
 # Get optimized solution variables
 x = prob.get_val("mesh.x_struct0", get_remote=True)[:-3:3]
 t_opt = prob["dv_struct"]

src/constitutive/TACSCompositeShellConstitutive.cpp

@@ -334,3 +334,22 @@ void TACSCompositeShellConstitutive::evalTangentHeatFlux(int elemIndex,
 const char *TACSCompositeShellConstitutive::getObjectName() {
   return constName;
 }
+
+/*
+  Get ply thicknesses
+*/
+void TACSCompositeShellConstitutive::getPlyThicknesses(
+    TacsScalar *_ply_thickness) {
+  for (int i = 0; i < num_plies; i++) {
+    _ply_thickness[i] = ply_thickness[i];
+  }
+}
+
+/*
+  Get ply angles
+*/
+void TACSCompositeShellConstitutive::getPlyAngles(TacsScalar *_ply_angles) {
+  for (int i = 0; i < num_plies; i++) {
+    _ply_angles[i] = ply_angles[i];
+  }
+}

src/constitutive/TACSCompositeShellConstitutive.h

@@ -75,6 +75,10 @@ class TACSCompositeShellConstitutive : public TACSShellConstitutive {
   // The name of the constitutive object
   const char *getObjectName();
 
+  // Get ply angles and thicknesses
+  void getPlyThicknesses(TacsScalar *_ply_thickness);
+  void getPlyAngles(TacsScalar *_ply_angles);
+
  private:
   // Store information about the design variable
   int num_plies;

src/elements/shell/TACSBeamElement.h

@@ -117,6 +117,12 @@ class TACSBeamRefAxisTransform : public TACSBeamTransform {
   }
   A2D::Vec3 &getRefAxis() { return axis; }
 
+  void getRefAxis(TacsScalar _axis[]) {
+    _axis[0] = axis.x[0];
+    _axis[1] = axis.x[1];
+    _axis[2] = axis.x[2];
+  }
+
  private:
   A2D::Vec3 axis;
   /* Tolerance for colinearity test in between beam axis and ref axis */

src/elements/shell/TACSShellElementTransform.h

@@ -107,6 +107,12 @@ class TACSShellRefAxisTransform : public TACSShellTransform {
     vec3Scale(invNorm, axis);
   }
 
+  void getRefAxis(TacsScalar _axis[]) {
+    _axis[0] = axis[0];
+    _axis[1] = axis[1];
+    _axis[2] = axis[2];
+  }
+
   void computeTransform(const TacsScalar Xxi[], const TacsScalar n0[],
                         TacsScalar T[]) {
     TacsScalar n[3];

src/elements/shell/TACSSpringElementTransform.h

@@ -47,6 +47,12 @@ class TACSSpringRefAxisTransform : public TACSSpringTransform {
     vec3Normalize(ref_dir);
   }
 
+  void getRefAxis(TacsScalar _axis[]) {
+    _axis[0] = ref_dir[0];
+    _axis[1] = ref_dir[1];
+    _axis[2] = ref_dir[2];
+  }
+
   /*
     Compute the local transformation from the global axis to the local
     reference frame for the spring element.
@@ -122,6 +128,16 @@ class TACSSpringRefFrameTransform : public TACSSpringTransform {
     crossProduct(&transform[6], &transform[0], &transform[3]);
   }
 
+  void getRefAxes(TacsScalar _axis_i[], TacsScalar _axis_j[]) {
+    _axis_i[0] = transform[0];
+    _axis_i[1] = transform[1];
+    _axis_i[2] = transform[2];
+
+    _axis_j[0] = transform[3];
+    _axis_j[1] = transform[4];
+    _axis_j[2] = transform[5];
+  }
+
   void computeTransform(const TacsScalar Xpts[], TacsScalar T[]) {
     memcpy(T, transform, 9 * sizeof(TacsScalar));
   }

tacs/TACS.pxd

@@ -100,6 +100,11 @@ cdef extern from "KSM.h":
         TACS_MAT_NORMAL
         TACS_MAT_TRANSPOSE
 
+cdef extern from "TACSMaterialProperties.h":
+    enum MaterialType:
+        TACS_ISOTROPIC_MATERIAL
+        TACS_ANISOTROPIC_MATERIAL
+
 # Special functions required for converting pointers
 cdef extern from "":
     TACSSchurMat* _dynamicSchurMat "dynamic_cast<TACSSchurMat*>"(TACSMat*)
@@ -357,6 +362,7 @@ cdef extern from "TACSElementModel.h":
 
 cdef class ElementModel:
     cdef TACSElementModel *ptr
+    cdef object con
 
 cdef inline _init_ElementModel(TACSElementModel *ptr):
     model = ElementModel()
@@ -387,6 +393,8 @@ cdef extern from "TACSElement.h":
 
 cdef class Element:
     cdef TACSElement *ptr
+    cdef object con
+    cdef object transform
 
 cdef inline _init_Element(TACSElement *ptr):
     elem = Element()
@@ -404,12 +412,18 @@ cdef class Function:
 cdef extern from "TACSConstitutive.h":
     cdef cppclass TACSConstitutive(TACSObject):
         int getNumStresses()
+        void evalStress(int, const double*, const TacsScalar*,
+                        const TacsScalar*, TacsScalar*)
+        TacsScalar evalDesignFieldValue(int, const double*,
+                             const TacsScalar*, int)
         void getFailureEnvelope(int, int, const double*,
                                 const TacsScalar*, const TacsScalar*,
                                 const TacsScalar*, TacsScalar*, TacsScalar*)
 
 cdef class Constitutive:
     cdef TACSConstitutive *ptr
+    cdef object props
+    cdef int nastranID
 
 cdef inline _init_Constitutive(TACSConstitutive *ptr):
     cons = Constitutive()

tacs/TACS.pyx

@@ -130,6 +130,10 @@ SEP_LARGEST = LARGEST
 SEP_SMALLEST_MAGNITUDE = SMALLEST_MAGNITUDE
 SEP_LARGEST_MAGNITUDE = LARGEST_MAGNITUDE
 
+# Import the material types
+ISOTROPIC_MATERIAL = TACS_ISOTROPIC_MATERIAL
+ANISOTROPIC_MATERIAL = TACS_ANISOTROPIC_MATERIAL
+
 # This wraps a C++ array with a numpy array for later useage
 cdef inplace_array_1d(int nptype, int dim1, void *data_ptr,
                       PyObject *ptr):
@@ -242,6 +246,7 @@ cdef class ElementBasis:
 cdef class ElementModel:
     def __cinit__(self, *args, **kwargs):
         self.ptr = NULL
+        self.con = None
         return
 
     def __dealloc__(self):
@@ -258,12 +263,19 @@ cdef class ElementModel:
             return self.ptr.getVarsPerNode()
         return 0
 
+    def getConstitutive(self):
+        if self.con:
+            return self.con
+        return None
+
 cdef class Element:
     """
     TACSElement base class
     """
     def __cinit__(self, *args, **kwargs):
         self.ptr = NULL
+        self.con = None
+        self.transform = None
         return
 
     def __dealloc__(self):
@@ -312,6 +324,16 @@ cdef class Element:
             return _init_ElementBasis(self.ptr.getElementBasis())
         return None
 
+    def getConstitutive(self):
+        if self.con:
+            return self.con
+        return None
+
+    def getTransform(self):
+        if self.transform:
+            return self.transform
+        return None
+
     def getElementType(self):
         if self.ptr:
             return self.ptr.getElementType()
@@ -449,6 +471,8 @@ cdef class Element:
 cdef class Constitutive:
     def __cinit__(self, *args, **kwargs):
         self.ptr = NULL
+        self.nastranID = 0
+        self.props = None
         return
 
     def __dealloc__(self):
@@ -467,6 +491,43 @@ cdef class Constitutive:
             return self.ptr.getNumStresses()
         return 0
 
+    def setNastranID(self, id):
+        """
+        Set property ID to be used in NASTRAN card for this object.
+        Should be set before `generateBDFCard` is called.
+
+        Args:
+            id (int): ID number to associate with this object's NASTRAN card
+        """
+        self.nastranID = id
+
+    def getNastranID(self):
+        """
+        Get property ID assigned in NASTRAN card for this object.
+
+        Returns:
+            id (int): ID number associated with this object's NASTRAN card
+        """
+        return self.nastranID
+
+    def getMaterialProperties(self):
+        """
+        Get the MaterialProperties class associated with this object
+
+        Returns:
+            prop (tacs.constitutive.MaterialProperties): TACS material property class associated with object.
+        """
+        return self.props
+
+    def generateBDFCard(self):
+        """
+        Generate pyNASTRAN card class based on current design variable values.
+
+        Returns:
+            card (pyNastran.bdf.cards.base_card.Property): pyNastran card holding property information
+        """
+        return None
+
     def getFailureEnvelope(self, sx, sy,
                            int elemIndex=0, int npts=100,
                            pt=[0.0, 0.0, 0.0], X=[0.0, 0.0, 0.0]):

tacs/constitutive.pxd

@@ -46,7 +46,9 @@ cdef extern from "TACSMaterialProperties.h":
                                TacsScalar, TacsScalar, TacsScalar,
                                TacsScalar, TacsScalar, TacsScalar)
         void setDensity(TacsScalar)
+        TacsScalar getDensity();
         void setSpecificHeat(TacsScalar)
+        void getIsotropicProperties(TacsScalar*, TacsScalar*)
         void getOrthotropicProperties(TacsScalar*, TacsScalar*, TacsScalar*,
                                       TacsScalar*, TacsScalar*, TacsScalar*,
                                       TacsScalar*, TacsScalar*, TacsScalar*)
@@ -56,6 +58,8 @@ cdef extern from "TACSMaterialProperties.h":
         void getCoefThermalExpansion(TacsScalar*, TacsScalar*, TacsScalar*)
         void getThermalConductivity(TacsScalar*, TacsScalar*, TacsScalar*)
 
+        MaterialType getMaterialType();
+
     cdef cppclass TACSOrthotropicPly(TACSObject):
         TACSOrthotropicPly(TacsScalar, TACSMaterialProperties*)
         void setKSWeight(TacsScalar)
@@ -64,6 +68,7 @@ cdef extern from "TACSMaterialProperties.h":
 
 cdef class MaterialProperties:
     cdef TACSMaterialProperties *ptr
+    cdef int nastranID
 
 cdef inline _init_MaterialProperties(TACSMaterialProperties *ptr):
     props = MaterialProperties()
@@ -115,6 +120,8 @@ cdef extern from "TACSCompositeShellConstitutive.h":
     cdef cppclass TACSCompositeShellConstitutive(TACSShellConstitutive):
         TACSCompositeShellConstitutive(int, TACSOrthotropicPly**, const TacsScalar*,
                                        const TacsScalar*, TacsScalar)
+        void getPlyThicknesses(TacsScalar*);
+        void getPlyAngles(TacsScalar*);
 
 cdef extern from "TACSLamParamShellConstitutive.h":
     cdef cppclass TACSLamParamShellConstitutive(TACSShellConstitutive):
@@ -148,6 +155,7 @@ cdef extern from "TACSIsoRectangleBeamConstitutive.h":
 cdef extern from "TACSGeneralMassConstitutive.h":
     cdef cppclass TACSGeneralMassConstitutive(TACSConstitutive):
         TACSGeneralMassConstitutive(const TacsScalar*)
+        void evalMassMatrix(int, const double*, const TacsScalar*, TacsScalar*)
 
 cdef class GeneralMassConstitutive(Constitutive):
     cdef TACSGeneralMassConstitutive *cptr
@@ -162,6 +170,7 @@ cdef extern from "TACSPointMassConstitutive.h":
 cdef extern from "TACSGeneralSpringConstitutive.h":
     cdef cppclass TACSGeneralSpringConstitutive(TACSConstitutive):
         TACSGeneralSpringConstitutive(TacsScalar*)
+        void evalMassMatrix(int, const double *, const TacsScalar *, TacsScalar *)
 
 cdef class GeneralSpringConstitutive(Constitutive):
     cdef TACSGeneralSpringConstitutive *cptr