Reorganize and add new examples

examples/fun3d_examples/README.md

@@ -1,38 +1,71 @@
-# FUN3D Examples #
-
-* `ate_wedge_optimization` - Optimization of a supersonic panel under supersonic flow using with aerothermoelastic analysis. Uses <font color="blue">FUN3D</font>, <font color="orange">TACS</font>, <font color="green">MELD</font>, `FuntofemNlbgs` driver.
-* `diamond_unsteady` - Unsteady forward and adjoint analysis of a diamond wedge wing structure under aeroelastic analysis. Uses <font color="blue">FUN3D</font>, <font color="orange">TACS</font>, <font color="green">MELD</font>.
-* `pyopt_togw_optimization` - Optimization of the CRM  aircraft structure under aeroelastic analysis. Uses <font color="blue">FUN3D</font>, <font color="orange">TACS</font>, <font color="green">MELD</font>, `FuntofemNlbgs` driver.
-* `sst_optimization` - Optimization of a supersonic transport wing under aerothermoelastic analysis using ksfailure and mass. Uses <font color="blue">FUN3D</font>, <font color="orange">TACS</font>, <font color="green">MELD</font>, `FuntofemNlbgs` driver.
-* `sst_unsteady` - Unsteady forward and adjoint analysis of a simplified supersonic transport wing geometry. Uses <font color="blue">FUN3D</font>, <font color="orange">TACS</font>, <font color="green">MELD</font>.
-
-### Supersonic Transport Wing ###
-<i>Directory</i> - `sst_optimization`
-The supersonic transport wing was the first demonstration of aerothermoelastic analysis with FUN3D and TACS on a realistic aircraft structure, included in the following paper. 
-```r
-Engelstad, S. P., Burke, B. J., Patel, R. N., Sahu, S., and Kennedy, G. J., “High-Fidelity Aerothermoelastic Optimization with
-Differentiable CAD Geometry,” AIAA Scitech 2023 Forum, National Harbor, MD, 2023. doi:10.2514/6.2023-0329.
+# FUN3D Examples
+
+Examples are organized by geometry group. Each group collects all FUN3D-based examples for a single aircraft or wing geometry.
+
+---
+
+## Super Simple Wing (SSW)
+
+The SSW is a simple rectangular wing geometry used for aeroelastic and aerothermal optimization studies.
+
+- **`ssw/aeroelastic_optimization/`** — Inviscid aeroelastic optimization of the Super Simple Wing. Uses FUN3D, TACS, and CAPS/ESP. Scripts 1–4 cover panel thickness optimization, shape optimization, and derivative testing.
+
+- **`ssw/ssw_meshdef_optimization/`** — Mesh deformation-based aeroelastic optimization of the Super Simple Wing. Uses FUN3D, TACS, and CAPS/ESP.
+
+---
+
+## Supersonic Transport Wing (SST)
+
+The SST is a realistic supersonic aircraft wing geometry used for aerothermoelastic optimization.
+
+The SST results images and README live directly in `sst/`. See `sst/sst_optimization/` for the complete runnable example.
+
+```
+Engelstad, S. P., Burke, B. J., Patel, R. N., Sahu, S., and Kennedy, G. J.,
+"High-Fidelity Aerothermoelastic Optimization with Differentiable CAD Geometry,"
+AIAA Scitech 2023 Forum, National Harbor, MD, 2023. doi:10.2514/6.2023-0329.
 ```
+
 <figure class="image">
-  <img src="sst_optimization/results/sst_opt_design.png" width=\linewidth/>
+  <img src="sst/results/sst_opt_design.png" width=600 />
   <figcaption><em>Optimal thicknesses for the supersonic transport wing.</em></figcaption>
 </figure>
 <br>
 <figure class="image">
-  <img src="sst_optimization/results/sst_fun3d_flow.png" width=600 />
+  <img src="sst/results/sst_fun3d_flow.png" width=600 />
   <figcaption><em>Pressure contours in the Mach 2.0 flow solved in FUN3D.</em></figcaption>
 </figure>
 
-### Computational Research Model ###
-<i>Directory</i> - `pyopt_togw_optimization`
-```r
-Jacobson, K., Kiviaho, J., Smith, M., and Kennedy, G., “An Aeroelastic Coupling Framework for Time-accurate Anal-
-ysis and Optimization,” 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2018.
-doi:10.2514/6.2018-0100.
+### `sst/sst_optimization/`
+
+Complete multi-step aerothermoelastic optimization example adapted from the research case. Includes mesh generation, sizing optimization, shape optimization, and fully coupled optimization scripts. Uses FUN3D, TACS, and CAPS/ESP.
+
 ```
-<figure class="image">
-  <img src="pyopt_togw_optimization/images/crm_thick_opt.png" width=400 />
-</figure> 
-<figure class="image">
-  <img src="pyopt_togw_optimization/images/crm_aero_struct_geom.png" width=400 />
-</figure>
+Engelstad, S. P., Burke, B. J., Patel, R. N., Sahu, S., and Kennedy, G. J.,
+"High-Fidelity Aerothermoelastic Optimization with Differentiable CAD Geometry,"
+AIAA Scitech 2023 Forum, National Harbor, MD, 2023. doi:10.2514/6.2023-0329.
+```
+
+---
+
+## Diamond
+
+- **`diamond/wedge_optimization/`** — Steady aerothermoelastic optimization of a supersonic diamond wedge panel. Minimizes average structural temperature subject to a mass constraint using a fully coupled FUN3D + TACS analysis. Uses a hand-built hexahedral BDF mesh (no CAPS/ESP required).
+
+---
+
+## Archive
+
+The `archive/` subdirectory preserves older examples that are no longer actively maintained. These examples are kept for historical reference but should not be used as templates for new work.
+
+- **`archive/ate_wedge_optimization/`** — Original version using deprecated APIs (`MassoudBody`, old `FUNtoFEMnlbgs` constructor, `PyOptOptimization`). Updated version is at `diamond/wedge_optimization/`.
+
+- **`archive/pyopt_togw_optimization/`** — Archived because it uses deprecated APIs (`MassoudBody`, `PyOptOptimization`). Original reference:
+
+  ```
+  Jacobson, K., Kiviaho, J., Smith, M., and Kennedy, G., "An Aeroelastic Coupling Framework for Time-accurate Analysis and Optimization," 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2018. doi:10.2514/6.2018-0100.
+  ```
+
+- **`archive/diamond_unsteady/`** — Archived because it is incomplete (no mesh files included).
+
+- **`archive/sst_unsteady/`** — Archived because it is incomplete (no mesh files included).

examples/fun3d_examples/diamond/wedge_optimization/README.md

@@ -0,0 +1,75 @@
+# Diamond Wedge — Aerothermoelastic Optimization
+
+Steady aerothermoelastic optimization of a supersonic diamond wedge panel. Minimizes average structural temperature subject to a mass constraint using a fully coupled FUN3D + TACS analysis.
+
+## Flow Conditions
+
+| Parameter | Value |
+|-----------|-------|
+| Mach | 0.5 |
+| V∞ | 148.67 m/s |
+| ρ∞ | 0.01037 kg/m³ |
+| q∞ | ~114.7 Pa |
+| T_ref | 300 K |
+| Wedge half-angle | 5° |
+
+## Design Variables
+
+10 panel thickness variables parameterized via a quadratic Bernstein basis (3 optimizer variables mapped to 10 TACS design variables in `gen_TACS_bdf_aero.py`).
+
+## Objective and Constraints
+
+- **Objective**: minimize average structural temperature
+- **Constraint**: structural mass ≤ 40 kg
+
+## Directory Layout
+
+```
+wedge_optimization/
+├── cfd/                    # FUN3D mesh and input files
+├── struct/                 # TACS BDF mesh (generated by gen_TACS_bdf_aero.py)
+├── design/                 # Optimization history and design variable files
+├── gen_TACS_bdf_aero.py    # Generates tacs_aero.bdf in struct/
+└── wedge_optimization.py   # Main optimization script
+```
+
+## Dependencies
+
+- FUN3D
+- TACS
+- mpi4py
+- pyoptsparse (SNOPT)
+
+## CFD Mesh
+
+No mesh is committed to this repository. The FUN3D mesh must be generated separately and placed in `cfd/`.
+
+**Geometry**: a 5° half-angle diamond wedge. The aerodynamic surface (FUN3D boundary tag 1) is the inclined upper face of the wedge:
+- Chord: x ∈ [1.0, 2.0] m (leading edge at x = 1 m, trailing edge at x = 2 m)
+- Span: y ∈ [−0.5, 0.5] m
+- Surface normal tilted 5° from vertical (matching the structural panel rotation in `gen_TACS_bdf_aero.py`)
+
+**Suggested mesh approach** using a structured mesh tool (e.g., Pointwise, ICEM CFD, or CGNS-based scripting):
+
+1. Define the wedge surface as a planar quad patch rotated 5° about the y-axis, with the leading edge at (1, ±0.5, 0) and trailing edge at (2, ±0.5, 0).
+2. Extrude a C-type or H-type block outward to a farfield boundary roughly 10–20 chord lengths away.
+3. Apply boundary conditions:
+   - `wall` (tag 1): the inclined wedge surface
+   - `farfield`: all remaining outer faces
+   - `symmetry` at y = ±0.5 m if exploiting span periodicity (optional)
+4. Export as a FUN3D-compatible unstructured grid (`.ugrid` or `.lb8.ugrid`) and place it in `cfd/` along with the FUN3D `input.nml`.
+
+For a quick inviscid mesh, a coarse H-block with ~50×20×20 cells is sufficient to demonstrate the coupled analysis. For a viscous case, add wall-normal clustering to achieve y⁺ ≈ 1 at the wall.
+
+## Run Order
+
+```bash
+# 1. Generate the structural mesh
+python gen_TACS_bdf_aero.py
+
+# 2. Generate the FUN3D CFD mesh (see CFD Mesh section above)
+#    and place the .ugrid and input.nml in cfd/
+
+# 3. Run the coupled optimization
+mpiexec -n <nprocs> python wedge_optimization.py
+```

examples/fun3d_examples/diamond/wedge_optimization/gen_TACS_bdf_aero.py

@@ -0,0 +1,151 @@
+"""
+gen_TACS_bdf_aero.py
+
+Generates the TACS BDF mesh for the diamond wedge aerothermoelastic example.
+Writes struct/tacs_aero.bdf relative to this script's directory.
+"""
+
+import numpy as np
+import os
+
+here = os.path.dirname(os.path.abspath(__file__))
+struct_dir = os.path.join(here, "struct")
+os.makedirs(struct_dir, exist_ok=True)
+
+nx = 11
+ny = 5
+nz = 5
+
+x_min = 1.0
+x_max = 2.0  # m
+y_min = -0.5
+y_max = 0.5  # m
+
+z_min = -0.015
+z_max = -0.005
+
+x = np.linspace(x_min, x_max, num=nx)
+y = np.linspace(y_min, y_max, num=ny)
+z = np.linspace(z_min, z_max, num=nz)
+theta = np.radians(5.0)
+nodes = np.arange(1, nx * ny * nz + 1, dtype=np.int64).reshape(nx, ny, nz)
+
+fp = open(os.path.join(struct_dir, "tacs_aero.bdf"), "w")
+fp.write("$ Input file for a rectangular plate\n")
+fp.write("SOL 103\nCEND\nBEGIN BULK\n")
+
+spclist = []
+spclistY = []
+spclistT = []
+# Write the grid points to a file
+for k in range(nz):
+    for j in range(ny):
+        for i in range(nx):
+            # Write the nodal data
+            spc = " "
+            coord_disp = 0
+            coord_id = 0
+            seid = 0
+
+            xpt = np.cos(theta) * x[i] - np.sin(theta) * z[k] + 1.0
+            ypt = y[j]
+            zpt = np.sin(theta) * x[i] + np.cos(theta) * z[k]
+
+            fp.write(
+                "%-8s%16d%16d%16.9e%16.9e*       \n"
+                % ("GRID*", nodes[i, j, k], coord_id, xpt, ypt)
+            )
+            fp.write(
+                "*       %16.9e%16d%16s%16d        \n" % (zpt, coord_disp, spc, seid)
+            )
+
+            # If the node is on one of the bottom outer edges of the plate,
+            # restrain it against displacement in any direction
+            if k == 0 and (i == 0 or i == nx - 1):
+                spclist.append(nodes[i, j, k])
+
+            # If the node is on one of the outer edges of the plate,
+            # restrain it against displacement in y direction
+            if k != 0 and (j == 0 or j == ny - 1):
+                spclistY.append(nodes[i, j, k])
+
+            # Set the temperature along the bottom edge of the plate
+            if k == 0:
+                spclistT.append(nodes[i, j, k])
+
+# Write out the linear hexahedral elements
+elem = 1
+for k in range(0, nodes.shape[2] - 1, 1):
+    for j in range(0, nodes.shape[1] - 1, 1):
+        for i in range(0, nodes.shape[0] - 1, 1):
+            # Set different part numbers for the elements on the
+            # lower and volume mesh
+            part_id = i + 1
+            if k == 0:
+                part_id = i + nodes.shape[0]
+            # Write the connectivity data
+            fp.write(
+                "%-8s%8d%8d%8d%8d%8d%8d%8d%8d\n%-8s%8d%8d\n"
+                % (
+                    "CHEXA",
+                    elem,
+                    part_id,
+                    nodes[i, j, k],
+                    nodes[i + 1, j, k],
+                    nodes[i + 1, j + 1, k],
+                    nodes[i, j + 1, k],
+                    nodes[i, j, k + 1],
+                    nodes[i + 1, j, k + 1],
+                    "*",
+                    nodes[i + 1, j + 1, k + 1],
+                    nodes[i, j + 1, k + 1],
+                )
+            )
+            elem += 1
+
+# # Add an extra layer of hexa elements and add the tractions through the
+# # TACSTraction3D class.
+# part_id = 2
+# k0 = 0
+# for j in range(0, nodes.shape[1]-1, 1):
+#     for i in range(0, nodes.shape[0]-1, 1):
+#         # Write the connectivity data
+#         fp.write('%-8s%8d%8d%8d%8d%8d%8d%8d%8d\n%-8s%8d%8d\n'%
+#                     ('CHEXA', elem, part_id,
+#                     nodes[i, j, k0], nodes[i+1, j ,k0],
+#                     nodes[i+1, j+1, k0], nodes[i, j+1, k0],
+#                     nodes[i, j, k0+1], nodes[i+1, j, k0+1], '*',
+#                     nodes[i+1, j+1, k0+1], nodes[i, j+1, k0+1]))
+#         elem += 1
+
+# Note: This approach will not work because the CQUAD4 pressure elements are not
+# yet implemented in TACS.
+#
+# # Write out the elements associated with the back-pressure. These are quads and
+# # are only associated with the back side of the structure
+# k0 = 0 # Set the plane of nodes where we will apply the surface traction
+# part_id = 2 # Set a different part id
+# for i in range(0, nodes.shape[0]-1, 1):
+#     for j in range(0, nodes.shape[1]-1, 1):
+#         n = [nodes[i, j, k0], nodes[i+1, j, k0],
+#              nodes[i+1, j+1, k0], nodes[i, j+1, k0]]
+#         # Note that the orientation of the element is reversed so that the normal
+#         # points in the -z direction
+#         fp.write('%-8s%8d%8d%8d%8d%8d%8d\n'%
+#                  ('CQUAD4', elem, part_id, n[0], n[3], n[2], n[1]))
+#         elem += 1
+
+for node in spclist:
+    spc = "123"
+    fp.write("%-8s%8d%8d%8s%8.6f\n" % ("SPC", 1, node, spc, 0.0))
+
+for node in spclistY:
+    spc = "2"
+    fp.write("%-8s%8d%8d%8s%8.6f\n" % ("SPC", 1, node, spc, 0.0))
+
+for node in spclistT:
+    spc = "4"
+    fp.write("%-8s%8d%8d%8s%8.4f\n" % ("SPC", 1, node, spc, 300.0))
+
+fp.write("END BULK")
+fp.close()