Deploying to gh-pages from @ ba94567 🚀

smdogroup · Oct 30, 2023 · 8292604 · 8292604
1 parent d5d42fb
commit 8292604
Show file tree

Hide file tree

Showing 16 changed files with 620 additions and 25 deletions.
diff --git a/_images/comp_opt_def.png b/_images/comp_opt_def.png
diff --git a/_images/pf_opt.png b/_images/pf_opt.png
diff --git a/_images/plate_pressure.png b/_images/plate_pressure.png
diff --git a/_sources/examples/Example-Beam_Optimization.rst.txt b/_sources/examples/Example-Beam_Optimization.rst.txt
@@ -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

diff --git a/_sources/examples/Example-Composite_Optimization.rst.txt b/_sources/examples/Example-Composite_Optimization.rst.txt
@@ -0,0 +1,239 @@
+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`.
+In this example a compliance minimization problem for a composite plate, 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 and the effects of stacking sequence. Instead, we homogenize or "smear" each ply angle's 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
+
+To begin we first 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 ``element_callback`` 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 ``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 ``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` class.
+
+.. 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)
+
+Fianlly, we define our :class:`~mphys.Multipoint` class.
+To do this, we instantiate the :class:`~tacs.mphys.builder.TacsBuilder` using the ``element_callback``, ``problem_setup``, and ``constraint_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")
+  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.
+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 optimization 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.
+
diff --git a/_sources/index.rst.txt b/_sources/index.rst.txt
@@ -26,6 +26,7 @@ Examples
    examples/Example-Plate
    examples/Example-Transient_Battery
    examples/Example-Beam_Optimization
+   examples/Example-Composite_Optimization
 
 References
 ==========

diff --git a/_sources/mphys/builder.rst.txt b/_sources/mphys/builder.rst.txt
@@ -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
 -------------

diff --git a/_sources/mphys/mphys.rst.txt b/_sources/mphys/mphys.rst.txt
@@ -2,9 +2,19 @@ MPhys
 *****
 `MPhys <https://pypi.org/project/mphys/>`_ is a package that standardizes high-fidelity multiphysics problems in OpenMDAO.
 MPhys provides a convenient interface for connecting TACS problems to optimizers and other physics disciplinary solvers for coupled analysis.
+MPhys is an optional TACS dependency, so in order to use the module within TACS it will need to be installed manually along with its dependencies.
+Assuming the user has already installed TACS, the only additional dependency the user will need to install would be `petsc4py <https://pypi.org/project/petsc4py/>`_.
+These packages can be installed easily in a conda environment using the commands below:
+
+::
+
+    pip install mphys
+    conda install -c conda-forge petsc4py
+
+
 The TACS MPhys interface consists of one main classes: a builder class called :class:`~tacs.mphys.builder.TacsBuilder`.
 For more information on the general MPhys interface users should see the MPhys `docs <https://openmdao.github.io/mphys/>`_.
-The details of the interfaces will be discussed in the sections below.
+The details of the :class:`~tacs.mphys.builder.TacsBuilder` interface will be discussed in the sections below.
 
 .. toctree::
   :maxdepth: 1

diff --git a/examples/Example-Beam_Optimization.html b/examples/Example-Beam_Optimization.html
@@ -15,7 +15,7 @@
     <script async="async" src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
     <link rel="index" title="Index" href="../genindex.html" />
     <link rel="search" title="Search" href="../search.html" />
-    <link rel="next" title="Theory" href="../theory/theory.html" />
+    <link rel="next" title="Composite plate optimization with MPhys" href="Example-Composite_Optimization.html" />
     <link rel="prev" title="Battery pack during thermal runaway" href="Example-Transient_Battery.html" /> 
   </head><body>
     <div class="related" role="navigation" aria-label="related navigation">
@@ -28,7 +28,7 @@ <h3>Navigation</h3>
           <a href="../py-modindex.html" title="Python Module Index"
              >modules</a> |</li>
         <li class="right" >
-          <a href="../theory/theory.html" title="Theory"
+          <a href="Example-Composite_Optimization.html" title="Composite plate optimization with MPhys"
              accesskey="N">next</a> |</li>
         <li class="right" >
           <a href="Example-Transient_Battery.html" title="Battery pack during thermal runaway"
@@ -109,7 +109,7 @@ <h1>Beam optimization with MPhys<a class="headerlink" href="#beam-optimization-w
     <span class="c1"># Set one thickness dv for every property group</span>
     <span class="n">con</span> <span class="o">=</span> <span class="n">constitutive</span><span class="o">.</span><span class="n">IsoRectangleBeamConstitutive</span><span class="p">(</span><span class="n">prop</span><span class="p">,</span> <span class="n">t</span><span class="o">=</span><span class="n">t</span><span class="p">,</span> <span class="n">w</span><span class="o">=</span><span class="n">w</span><span class="p">,</span> <span class="n">tNum</span><span class="o">=</span><span class="n">dvNum</span><span class="p">)</span>
 
-    <span class="c1"># Defines local y/width direction for beam</span>
+    <span class="c1"># Defines local y/thickness direction for beam</span>
     <span class="n">refAxis</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">array</span><span class="p">([</span><span class="mf">0.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">,</span> <span class="mf">0.0</span><span class="p">])</span>
     <span class="n">transform</span> <span class="o">=</span> <span class="n">elements</span><span class="o">.</span><span class="n">BeamRefAxisTransform</span><span class="p">(</span><span class="n">refAxis</span><span class="p">)</span>
 
@@ -226,12 +226,16 @@ <h1>Beam optimization with MPhys<a class="headerlink" href="#beam-optimization-w
 <span class="c1"># Get analytical solution</span>
 <span class="n">t_exact</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">6</span> <span class="o">*</span> <span class="p">(</span><span class="n">L</span> <span class="o">-</span> <span class="n">x</span><span class="p">)</span> <span class="o">*</span> <span class="n">V</span> <span class="o">/</span> <span class="n">w</span> <span class="o">/</span> <span class="n">ys</span><span class="p">)</span>
 
+<span class="c1"># Compute max thickness value</span>
+<span class="n">t0</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mi">6</span> <span class="o">*</span> <span class="n">L</span> <span class="o">*</span> <span class="n">V</span> <span class="o">/</span> <span class="n">w</span> <span class="o">/</span> <span class="n">ys</span><span class="p">)</span>
+
 <span class="c1"># Plot results for solution</span>
-<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">t_opt</span><span class="p">,</span> <span class="s2">&quot;o&quot;</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">t_exact</span><span class="p">)</span>
+<span class="n">plt</span><span class="o">.</span><span class="n">plot</span><span class="p">(</span><span class="n">x</span> <span class="o">/</span> <span class="n">L</span><span class="p">,</span> <span class="n">t_opt</span> <span class="o">/</span> <span class="n">t0</span><span class="p">,</span> <span class="s2">&quot;o&quot;</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">t_exact</span> <span class="o">/</span> <span class="n">t0</span><span class="p">)</span>
 <span class="n">plt</span><span class="o">.</span><span class="n">legend</span><span class="p">([</span><span class="s2">&quot;optimized&quot;</span><span class="p">,</span> <span class="s2">&quot;analytical&quot;</span><span class="p">])</span>
-<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="s2">&quot;t(x)&quot;</span><span class="p">)</span>
-<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="s2">&quot;x&quot;</span><span class="p">)</span>
+<span class="n">plt</span><span class="o">.</span><span class="n">ylabel</span><span class="p">(</span><span class="sa">r</span><span class="s2">&quot;$\frac{t(x)}</span><span class="si">{t_0}</span><span class="s2">$&quot;</span><span class="p">,</span> <span class="n">fontsize</span><span class="o">=</span><span class="mi">16</span><span class="p">)</span>
+<span class="n">plt</span><span class="o">.</span><span class="n">xlabel</span><span class="p">(</span><span class="sa">r</span><span class="s2">&quot;$\frac</span><span class="si">{x}{L}</span><span class="s2">$&quot;</span><span class="p">,</span> <span class="n">fontsize</span><span class="o">=</span><span class="mi">16</span><span class="p">,</span> <span class="n">labelpad</span><span class="o">=-</span><span class="mi">5</span><span class="p">)</span>
 <span class="n">plt</span><span class="o">.</span><span class="n">title</span><span class="p">(</span><span class="s2">&quot;Optimal beam thickness profile&quot;</span><span class="p">)</span>
+<span class="n">plt</span><span class="o">.</span><span class="n">text</span><span class="p">(</span><span class="mf">0.05</span><span class="p">,</span> <span class="mf">0.25</span><span class="p">,</span> <span class="sa">r</span><span class="s2">&quot;$t_0 = \sqrt{\frac</span><span class="si">{6VL}</span><span class="s2">{w\sigma_y}}$&quot;</span><span class="p">,</span> <span class="n">fontsize</span><span class="o">=</span><span class="mi">12</span><span class="p">)</span>
 <span class="n">plt</span><span class="o">.</span><span class="n">show</span><span class="p">()</span>
 </pre></div>
 </div>
@@ -252,8 +256,8 @@ <h4>Previous topic</h4>
   </div>
   <div>
     <h4>Next topic</h4>
-    <p class="topless"><a href="../theory/theory.html"
-                          title="next chapter">Theory</a></p>
+    <p class="topless"><a href="Example-Composite_Optimization.html"
+                          title="next chapter">Composite plate optimization with MPhys</a></p>
   </div>
 <div id="searchbox" style="display: none" role="search">
   <h3 id="searchlabel">Quick search</h3>
@@ -279,7 +283,7 @@ <h3>Navigation</h3>
           <a href="../py-modindex.html" title="Python Module Index"
              >modules</a> |</li>
         <li class="right" >
-          <a href="../theory/theory.html" title="Theory"
+          <a href="Example-Composite_Optimization.html" title="Composite plate optimization with MPhys"
              >next</a> |</li>
         <li class="right" >
           <a href="Example-Transient_Battery.html" title="Battery pack during thermal runaway"