Mphys helper functions (#323)

* Add mphys helper function for adding constraints

* Alter composite optimisation example to use new mphys helper function

* Add a function decorator for logging function calls to TACS MPhys components

* Revert "Add a function decorator for logging function calls to TACS MPhys components"

This reverts commit e4e9363.

* Remove unnecessary model input

* Test `add_tacs_constraints` in mphys integration test

* Need to add a nonlinear constraint to get around OpenMDAO bug

* Add an objective too for realism

examples/plate/optimize_composite_plate.py

@@ -1,6 +1,7 @@
 """
 Mass minimization of uCRM wingbox subject to a constant vertical force
 """
+
 import os
 
 import openmdao.api as om
@@ -10,6 +11,7 @@
 
 from tacs import elements, constitutive, functions
 from tacs.mphys import TacsBuilder
+from tacs.mphys.utils import add_tacs_constraints
 
 # BDF file containing mesh
 bdf_file = os.path.join(os.path.dirname(__file__), "partitioned_plate.bdf")
@@ -42,7 +44,9 @@
 
 
 # Callback function used to setup TACS element objects and DVs
-def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs):
+def element_callback(
+    dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs
+):
     # Create ply object
     ortho_prop = constitutive.MaterialProperties(
         rho=rho,
@@ -105,7 +109,12 @@ def constraint_setup(scenario_name, fea_assembler, constraint_list):
     constr = fea_assembler.createDVConstraint("ply_fractions")
     allComponents = fea_assembler.selectCompIDs()
     constr.addConstraint(
-        "sum", allComponents, dvIndices=[0, 1, 2, 3], dvWeights=[1.0, 1.0, 1.0, 1.0]
+        "sum",
+        allComponents,
+        dvIndices=[0, 1, 2, 3],
+        dvWeights=[1.0, 1.0, 1.0, 1.0],
+        lower=1.0,
+        upper=1.0,
     )
     constraint_list.append(constr)
 
@@ -126,14 +135,21 @@ def setup(self):
         dvs = self.add_subsystem("dvs", om.IndepVarComp(), promotes=["*"])
         dvs.add_output("dv_struct", dv_array)
 
-        self.add_subsystem("mesh", struct_builder.get_mesh_coordinate_subsystem())
+        self.add_subsystem(
+            "mesh", struct_builder.get_mesh_coordinate_subsystem()
+        )
         self.mphys_add_scenario(
             "pressure_load", ScenarioStructural(struct_builder=struct_builder)
         )
-        self.mphys_connect_scenario_coordinate_source("mesh", "pressure_load", "struct")
+        self.mphys_connect_scenario_coordinate_source(
+            "mesh", "pressure_load", "struct"
+        )
 
         self.connect("dv_struct", "pressure_load.dv_struct")
 
+    def configure(self):
+        add_tacs_constraints(self.pressure_load)
+
 
 ################################################################################
 # OpenMDAO setup
@@ -146,10 +162,11 @@ def setup(self):
 # Declare design variables, objective, and constraint
 model.add_design_var("dv_struct", lower=0.0, upper=1.0)
 model.add_objective("pressure_load.compliance")
-model.add_constraint("pressure_load.ply_fractions.sum", equals=1.0, linear=True)
 
 # Configure optimizer
-prob.driver = om.ScipyOptimizeDriver(debug_print=["objs", "nl_cons"], maxiter=100)
+prob.driver = om.ScipyOptimizeDriver(
+    debug_print=["objs", "nl_cons", "ln_cons"], maxiter=100
+)
 prob.driver.options["optimizer"] = "SLSQP"
 
 # Setup OpenMDAO problem

tacs/mphys/__init__.py

@@ -1,3 +1,4 @@
 from tacs.mphys.builder import TacsBuilder
+import tacs.mphys.utils
 
-__all__ = ["TacsBuilder"]
+__all__ = ["TacsBuilder", "utils"]

tacs/mphys/utils.py

@@ -0,0 +1,62 @@
+import warnings
+import tacs
+
+
+def add_tacs_constraints(scenario):
+    """Call this in the configure method to add all TACS constraints functions
+    in a TacsPostcouplingGroup as constraints in an OpenMDAO model. This saves
+    you having to call OpenMDAO's `add_constraint` method for each constraint
+    function yourself.
+
+    This function relies on you having defined the lower and/or upper bounds
+    of the constraints when you created them in your `constraint_setup`
+    function.
+
+    Parameters
+    ----------
+    scenario : MPhys scenario
+        Scenario containing a TACS postcoupling group with constraints
+    """
+    # First, verify that the model has a struct_post component that is a
+    # TacsPostcouplingGroup
+    if not hasattr(scenario, "struct_post"):
+        raise ValueError(
+            f"Scenario {scenario.name} does not have a struct_post component"
+        )
+    if not isinstance(
+        scenario.struct_post, tacs.mphys.postcoupling.TacsPostcouplingGroup
+    ):
+        raise ValueError(
+            f"{scenario.name}.struct_post is not a TacsPostcouplingGroup, cannot add TACS constraints"
+        )
+    # If the provided scenario has no constraints component, warn the user and return
+    if not hasattr(scenario.struct_post, "constraints"):
+        warnings.warn(
+            f"Scenario {scenario.name} contains no TACS constraints to add",
+            stacklevel=2,
+        )
+        return
+    else:
+        for system in scenario.struct_post.constraints.system_iter():
+            constraint = system.constr
+            constraintFuncNames = constraint.getConstraintKeys()
+            bounds = {}
+            constraint.getConstraintBounds(bounds)
+            for conName in constraintFuncNames:
+                if scenario.comm.rank == 0:
+                    print("Adding TACS constraint: ", conName)
+                name = f"{system.name}.{conName}"
+                # Panel length constraints are nonlinear, all other constrain
+                # types (DV and adjacency) are linear
+                if isinstance(
+                    constraint,
+                    tacs.constraints.panel_length.PanelLengthConstraint,
+                ):
+                    scenario.add_constraint(name, equals=0.0, scaler=1.0)
+                else:
+                    lb = bounds[f"{system.name}_{conName}"][0]
+                    ub = bounds[f"{system.name}_{conName}"][1]
+                    if all(lb == ub):
+                        scenario.add_constraint(name, equals=lb, linear=True)
+                    else:
+                        scenario.add_constraint(name, lower=lb, upper=ub, linear=True)

tests/integration_tests/test_mphys_struct_plate.py

@@ -33,6 +33,10 @@
 # KS function weight
 ksweight = 10.0
 
+# Adjacency constraint bounds
+adj_lb = -2.5e-3
+adj_ub = 2.5e-3
+
 
 class ProblemTest(OpenMDAOTestCase.OpenMDAOTest):
     N_PROCS = 2  # this is how many MPI processes to use for this TestCase.
@@ -105,7 +109,7 @@ def constraint_setup(scenario_name, fea_assembler, constraints):
             """
             # Setup adjacency constraints for panel thicknesses
             constr = fea_assembler.createAdjacencyConstraint("adjacency")
-            constr.addConstraint("PANEL")
+            constr.addConstraint("PANEL", lower=-2.5e-3, upper=2.5e-3)
             constr_list = [constr]
             return constr_list
 
@@ -145,6 +149,11 @@ def setup(self):
                 self.connect("dv_struct", "analysis.dv_struct")
                 self.connect("f_struct", "analysis.f_struct")
 
+            def configure(self):
+                tacs.mphys.utils.add_tacs_constraints(self.analysis)
+                self.add_constraint("analysis.ks_vmfailure", upper=1.0)
+                self.add_objective("analysis.mass")
+
         prob = om.Problem()
         prob.model = Top()
 
@@ -157,12 +166,23 @@ def setup_funcs(self):
         """
         return FUNC_REFS, wrt
 
+    # def test_add_tacs_constraints(self):
+    #     # prob = self.setup_problem(dtype=float)
+    #     # prob.setup()
+    #     prob = self.prob
+    #     constraints = prob.model.get_constraints()
+    #     self.assertIn("analysis.adjacency.PANEL", constraints)
+    #     adjacency = constraints["analysis.adjacency.PANEL"]
+    #     self.assertTrue(adjacency["linear"])
+    #     lower_bound = adjacency["lower"]
+    #     upper_bound = adjacency["upper"]
+    #     np.testing.assert_equal(lower_bound, adj_lb)
+    #     np.testing.assert_equal(upper_bound, adj_ub)
+
     def test_get_tagged_indices(self):
         """
         Test the get_tagged_indices method
         """
-        prob = self.setup_problem(dtype=float)
-        prob.setup()
 
         # We want to test that:
         # - For each comp_id, get_tagged_indices returns the same indices as `getLocalNodeIDsForComps`
@@ -192,3 +212,9 @@ def test_get_tagged_indices(self):
         trueNodeIDs = FEAAssembler.getLocalNodeIDsForComps([compIDs[0], compIDs[-1]])
         taggedIndIDs = self.tacs_builder.get_tagged_indices(tags)
         self.assertEqual(sorted(trueNodeIDs), sorted(taggedIndIDs))
+
+
+if __name__ == "__main__":
+    import unittest
+
+    unittest.main()