Adding composite opt example

docs/source/examples/Example-Beam_Optimization.rst

@@ -76,7 +76,7 @@ assign a design variable number for the thickness parameter, and return a :class
       # Set one thickness dv for every property group
       con = constitutive.IsoRectangleBeamConstitutive(prop, t=t, w=w, tNum=dvNum)
 
-      # Defines local y/width direction for beam
+      # Defines local y/thickness direction for beam
       refAxis = np.array([0.0, 1.0, 0.0])
       transform = elements.BeamRefAxisTransform(refAxis)
 
@@ -203,12 +203,16 @@ Finally, we can plot the optimized thickness distribution using matplotlib and c
   # Get analytical solution
   t_exact = np.sqrt(6 * (L - x) * V / w / ys)
 
+  # Compute max thickness value
+  t0 = np.sqrt(6 * L * V / w / ys)
+
   # Plot results for solution
-  plt.plot(x, t_opt, "o", x, t_exact)
+  plt.plot(x / L, t_opt / t0, "o", x, t_exact / t0)
   plt.legend(["optimized", "analytical"])
-  plt.ylabel("t(x)")
-  plt.xlabel("x")
+  plt.ylabel(r"$\frac{t(x)}{t_0}$", fontsize=16)
+  plt.xlabel(r"$\frac{x}{L}$", fontsize=16, labelpad=-5)
   plt.title("Optimal beam thickness profile")
+  plt.text(0.05, 0.25, r"$t_0 = \sqrt{\frac{6VL}{w\sigma_y}}$", fontsize=12)
   plt.show()
 
 .. image:: images/beam_plot.png

docs/source/examples/Example-Composite_Optimization.rst

@@ -0,0 +1,245 @@
+Composite plate optimization with MPhys
+***************************************
+.. note:: The script for this example can be found under the `examples/plate/` directory.
+
+This example further demonstrates TACS structural optimization capabilities through :ref:`mphys/mphys:MPhys`.
+The beam model that we will be using for this problem is a rectangular beam,
+cantilevered, with a shear load applied at the tip. The beam is discretized using
+100 beam elements along it's span.
+
+In this example a composite plate compliance minimization problem, based off of one originally proposed by Lund and Stegmann [`1`_], is solved.
+A 1 m x 1 m x 0.05 m composite plate is clamped on all edges and subjected to
+a uniform pressure of 100 kPa loading on the top.
+The plate is discretized into 100 quad shell elements each having its own independent laminate layup.
+The goal of this example will be to optimize the layup of each element in the plate in order to minimize the total compliance energy.
+A diagram of the problem is shown below.
+
+.. image:: images/plate_pressure.png
+  :width: 800
+  :alt: Plate problem
+
+To make this problem tractable for gradient-based optimization we will make the following simplifications:
+
+  1. For each layup we can ony select plies of the following angles: :math:`0^\circ, 45^\circ, -45^\circ, 90^\circ`.
+
+  2. We will neglect the discrete nature of the plies the effects of stacking sequence and instead homogenize or "smear" each ply angles contribution to laminate stiffness based on its proportion of plies in the layup or "ply fraction".
+
+The optimization problem can now be summarized as follows:
+
+.. image:: images/comp_opt_def.png
+  :width: 800
+  :alt: Optimization problem
+
+First, we import required libraries, define the model bdf file, and define important problem constants:
+
+.. code-block:: python
+
+  import os
+
+  import openmdao.api as om
+  import numpy as np
+  from mphys import Multipoint
+  from mphys.scenario_structural import ScenarioStructural
+
+  from tacs import elements, constitutive, functions
+  from tacs.mphys import TacsBuilder
+
+  # BDF file containing mesh
+  bdf_file = os.path.join(os.path.dirname(__file__), "partitioned_plate.bdf")
+
+  # Material properties
+  rho = 1550.0
+  E1 = 54e9
+  E2 = 18e9
+  nu12 = 0.25
+  G12 = 9e9
+  G13 = 9e9
+  Xt = 2410.0e6
+  Xc = 1040.0e6
+  Yt = 73.0e6
+  Yc = 173.0e6
+  S12 = 71.0e6
+
+  # Shell thickness
+  ply_thickness = 1.25e-3  # m
+  plate_thickness = 0.05  # m
+  tMin = 0.002  # m
+  tMax = 0.05  # m
+
+  # Ply angles/initial ply fractions
+  ply_angles = np.deg2rad([0.0, 45.0, -45.0, 90.0])
+  ply_fractions = np.array([0.25, 0.25, 0.25, 0.25])
+
+  # Pressure load to apply to plate
+  P = 100e3
+
+Next we define an :func:`~tacs.pytacs.elemCallBack` function for setting up the TACS elements and design variables.
+We use the :class:`~tacs.constitutive.SmearedCompositeShellConstitutive` class here for the constitutive properties, and
+assign four design variable numbers to each element (one for each ply fraction), and return a :class:`~tacs.elements.Quad4Shell` element class.
+
+.. code-block:: python
+
+  # Callback function used to setup TACS element objects and DVs
+  def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs):
+      # Create ply object
+      ortho_prop = constitutive.MaterialProperties(
+          rho=rho,
+          E1=E1,
+          E2=E2,
+          nu12=nu12,
+          G12=G12,
+          G13=G13,
+          G23=G13,
+          Xt=Xt,
+          Xc=Xc,
+          Yt=Yt,
+          Yc=Yc,
+          S12=S12,
+      )
+      ortho_ply = constitutive.OrthotropicPly(ply_thickness, ortho_prop)
+      # Create the layup list (one for each angle)
+      ortho_layup = [ortho_ply, ortho_ply, ortho_ply, ortho_ply]
+      # Assign each ply fraction a unique DV
+      ply_fraction_dv_nums = np.array(
+          [dvNum, dvNum + 1, dvNum + 2, dvNum + 3], dtype=np.intc
+      )
+      # Create smeared stiffness object based on ply angles/fractions
+      con = constitutive.SmearedCompositeShellConstitutive(
+          ortho_layup,
+          plate_thickness,
+          ply_angles,
+          ply_fractions,
+          ply_fraction_dv_nums=ply_fraction_dv_nums,
+      )
+
+      # Define reference axis to define local 0 deg direction
+      refAxis = np.array([1.0, 0.0, 0.0])
+      transform = elements.ShellRefAxisTransform(refAxis)
+
+      # Pass back the appropriate tacs element object
+      elem = elements.Quad4Shell(transform, con)
+
+      return elem
+
+We define a :func:`problem_setup` to add fixed loads and eval functions.
+Here we specify the plate compliance energy (:class:`~tacs.functions.Compliance`) as an output for our analysis
+and add our 100 kPa pressure load.
+
+.. code-block:: python
+
+  def problem_setup(scenario_name, fea_assembler, problem):
+      """
+      Helper function to add fixed forces and eval functions
+      to structural problems used in tacs builder
+      """
+
+      # Add TACS Functions
+      problem.addFunction("compliance", functions.Compliance)
+
+      # Add forces to static problem
+      allComponents = fea_assembler.selectCompIDs()
+      problem.addPressureToComponents(allComponents, P)
+
+For our last helper function we define a :func:`constraint_setup` function.
+This function can be used to add additional relational constraints to the design variables we defined in the ``element_callback``.
+In particular, we want to enforce a new constraint (100 in total) such that the ply fractions within each element should sum to unity.
+We can accomplish this by utilizing the :class:`~tacs.constraints.DVConstraint`.
+
+.. code-block:: python
+
+  def constraint_setup(scenario_name, fea_assembler, constraint_list):
+      """
+      Helper function to setup tacs constraint classes
+      """
+      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]
+      )
+      constraint_list.append(constr)
+
+Here we define our :class:`~mphys.Multipoint` (essentially an OpenMDAO ``Group``) which will contain our analysis :class:`~mphys.Scenario`.
+To do this, we instantiate the :class:`~tacs.mphys.builder.TacsBuilder` using the ``element_callback`` and ``problem_setup`` we defined above.
+We create OpenMDAO ``Component``'s to feed design variable and mesh inputs to the ``Scenario`` component.
+We use this builder to create an MPhys :class:`~mphys.StructuralScenario`.
+
+.. code-block:: python
+
+  class PlateModel(Multipoint):
+      def setup(self):
+          struct_builder = TacsBuilder(
+              mesh_file=bdf_file,
+              element_callback=element_callback,
+              problem_setup=problem_setup,
+              constraint_setup=constraint_setup,
+              coupled=False,
+              check_partials=True,
+          )
+          struct_builder.initialize(self.comm)
+          dv_array = struct_builder.get_initial_dvs()
+
+          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.mphys_add_scenario(
+              "pressure_load", ScenarioStructural(struct_builder=struct_builder)
+          )
+          self.mphys_connect_scenario_coordinate_source("mesh", "pressure_load", "struct")
+
+          self.connect("dv_struct", "pressure_load.dv_struct")
+
+At this point we setup the OpenMDAO ``Problem`` class that we will use to perform our optimization.
+We assign our ``PlateModel`` to the problem class and set ``ScipyOptimizeDriver``.
+We define our design variables, constraint, and objective.
+Finally we run the problem driver to optimize the problem.
+
+.. code-block:: python
+
+  prob = om.Problem()
+  prob.model = PlateModel()
+  model = prob.model
+
+  # Declare design variables, objective, and constraint
+  model.add_design_var("dv_struct", lower=0.0, upper=1.0)
+  model.add_objective("pressure_load.compliance", scaler=1.0)
+  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.options["optimizer"] = "SLSQP"
+
+  # Setup OpenMDAO problem
+  prob.setup()
+
+  # Output N2 representation of OpenMDAO model
+  om.n2(prob, show_browser=False, outfile="tacs_struct.html")
+
+  # Run optimization
+  prob.run_driver()
+
+After the optimization completes the user should see a print out to screen like shown below.
+
+>>> Optimization terminated successfully    (Exit mode 0)
+>>>             Current function value: 8.571649588963465
+>>>             Iterations: 34
+>>>             Function evaluations: 34
+>>>             Gradient evaluations: 34
+>>> Optimization Complete
+>>> -----------------------------------
+
+Once the optimization is complete we can post-process results.
+We can write our optimized beam model to a BDF file so they can
+be processed in other commonly used FEM software.
+The ``f5`` solution file at each optimization iteration can also be converted to a Tecplot or Paraview files using ``f5totec`` or ``f5tovtk``, respectively.
+The optimized ply fraction distributions for each angle can be visualized by plotting the contours of the following variables in Tecplot or Paraview: ``dv2``, ``dv3``, ``dv4``, ``dv5``.
+A visualization of the optimized result is shown below:
+
+.. image:: images/pf_opt.png
+  :width: 800
+  :alt: Plate solution
+
+.. rubric:: References
+
+.. [1] Lund, E. and Stegmann, J., “On structural optimization of composite shell structures using a discrete constitutive parametrization,” Wind Energy, Vol. 8, No. 1, 2005, pp. 109–124.
+

docs/source/index.rst

@@ -26,6 +26,7 @@ Examples
    examples/Example-Plate
    examples/Example-Transient_Battery
    examples/Example-Beam_Optimization
+   examples/Example-Composite_Optimization
 
 References
 ==========

docs/source/mphys/builder.rst

@@ -3,7 +3,6 @@ TacsBuilder class
 TACS MPhys interface centers around the :class:`~tacs.mphys.builder.TacsBuilder` class.
 The :class:`~tacs.mphys.builder.TacsBuilder` is responsible for creating :ref:`Scenarios <mphys:scenario_groups>`,
 an OpenMDAO group that contains an analysis condition in a multipoint optimization.
-:ref:`mphys:builders` .
 
 API Reference
 -------------

examples/beam/beam_opt.py

@@ -56,7 +56,7 @@ def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **k
     # Set one thickness dv for every property group
     con = constitutive.IsoRectangleBeamConstitutive(prop, t=t, w=w, tNum=dvNum)
 
-    # Defines local y/width direction for beam
+    # Defines local y/thickness direction for beam
     refAxis = np.array([0.0, 1.0, 0.0])
     transform = elements.BeamRefAxisTransform(refAxis)