Implemented multi-section surface parameterization (#446)

* Feature: Start port of geomEngine3 from Matlab to OSA.WIP

* Implemented mesh generator for multisection wings

* Added test code for plotting mesh

* Fixed bugs related to chordwise panelling not working properly at wing  root. Started proper  multi-section surface implementation.

* Complete support for multi-section wings

* Clean-up plotting functions

* Removed stray files

* Fixed plotting bugs

* Completed mesh plotting code

* Add OAS mesh output

* Start separating code into helper functions

* Continue work on helper functions.

* Code clean-up

* Main mesh generator functions implemented and commented

* Complete the initial implementation of the user custom mesh utility

* Code cleanup

* Fix flake8 and black versions used

* Fix remaining formatting issues

* Start new Geometry group implementation

* Add support for sectional mesh generation

* Enabled basic section adding functionality

* Implemnented multi section geom group

* First set of mesh bug fixes

* Fixed implementation of mesh generation and multi section mesh geometry group.

* Implemented basic continuity constraint

* Implemented support for Aero Point on multisection. Requires a vortex mesh bug fix to function properly with symmetry.

* Introduced the multi-section geometry group and associated components.

* Refactored code, asymmetry support, and C0 by constuction support

* Clean-up and added comments

* Fixed t_over_c b-splint for multi section

* Added tests in temp location

* Cleaned up section dictionary handling

* Improved dictionary handling. Added support for multi-surface + multi-section

* Added multi-section geometry unit tests

* Regression tests implemented

* aero group updates

* Remove rouge print statements and remove accidental use of finite difference partials for mesh unification

* Added docs

* remove personal test files

* Corrected tests location

* Formatting updates and reorganization of tests

* Fix one more formatting issue

* Added single and multisection comparison test

* More black and flake8 fixes

* Should fix the black and flake8

* Flake8 fixes

* Fixed doc navigation

* Formatting changes for change request

* Update test initial design

* Reworked the constraint vs construction test

* Cleaned up leftover stuff

* Fix flake8

* Disable SNOPT test for now

* Fix failing test

* SNOPT ON again

* Go back to SLSQP

* Reorganized test to help it pass

* Try to fix test again

* Fix formatting

* Removed SNOPT test

* Update test util names

* Start moving test surfaces to testing utils functions

* Reorganized and reworked all tests

* Update multi-section docs and examples

* Flake8 fixes

* Merge check surface keys into one function

* reformat function and variable names for PEP8

* Update surface keys user reference

* Cleaned up test in response to feedback

* Fix unit tests

* Bump minor version

openaerostruct/__init__.py

@@ -1 +1 @@
-__version__ = "2.8.0"
+__version__ = "2.9.0"

openaerostruct/aerodynamics/aero_groups.py

@@ -30,6 +30,37 @@ def setup(self):
         surfaces = self.options["surfaces"]
         rotational = self.options["rotational"]
 
+        # Check for multi-section surfaces and create suitable surface dictionaries for them
+        for i, surface in enumerate(surfaces):
+            # If multisection mesh then build a single surface with the unified mesh data
+            if "is_multi_section" in surface.keys():
+                import copy
+
+                target_keys = [
+                    # Essential Info
+                    "name",
+                    "symmetry",
+                    "S_ref_type",
+                    "ref_axis_pos",
+                    "mesh",
+                    # aerodynamics
+                    "CL0",
+                    "CD0",
+                    "with_viscous",
+                    "with_wave",
+                    "groundplane",
+                    "k_lam",
+                    "t_over_c_cp",
+                    "c_max_t",
+                ]
+
+                # Constructs a surface dictionary and adds the specified supported keys and values from the mult-section surface dictionary.
+                aeroSurface = {}
+                for k in set(surface).intersection(target_keys):
+                    aeroSurface[k] = surface[k]
+                # print(aeroSurface["name"])
+                surfaces[i] = copy.deepcopy(aeroSurface)
+
         # Loop through each surface and connect relevant parameters
         for surface in surfaces:
             name = surface["name"]

openaerostruct/docs/advanced_features.rst

@@ -12,6 +12,7 @@ These examples show some advanced features in OpenAeroStruct.
    advanced_features/geometry_manipulation.rst
    advanced_features/custom_mesh_example.rst
    advanced_features/openvsp_mesh_example.rst
+   advanced_features/multi_section_surfaces.rst
    advanced_features/multiple_surfaces.rst
    advanced_features/multipoint.rst
    advanced_features/multipoint_parallel.rst

openaerostruct/docs/advanced_features/multi_section_surfaces.rst

@@ -0,0 +1,107 @@
+.. _Custom_Mesh:
+
+Multi-section Surfaces
+==========================
+
+OpenAeroStruct features the ability to specify surfaces as a series of sequentially connected sections. 
+Rather than controling the geometry of the surface as a whole, the optimizer can control the geometric parameters of each section individually using any of the geometric transformations available in OpenAeroStruct. 
+This feature was developed for the modular wing morphing applications but can be useful in situations where the user wishes to optimize a particular wing section while leaving the others fixed. 
+
+*Please note that aerostructural optimization with multi-section wings is currently not supported*
+
+This example script demonstrate the usage of the multi-section wing geometry features in OpenAeroStruct.
+We first start with the induced drag minimization of a simple two-section symmetrical wing.
+
+Let's start with the necessary imports.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec.py
+   :start-after: checkpoint 0
+   :end-before: checkpoint 1
+
+
+Then we are ready to discuss the multi-section parameterization and multi-section surface dictionary.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec.py
+   :start-after: checkpoint 1
+   :end-before: checkpoint 2
+
+Next, we setup the flow variables.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec.py
+   :start-after: checkpoint 2
+   :end-before: checkpoint 3
+
+
+Giving the optimizer control of each wing section without any measure to constrain its motion can lead to undesirable geometries, such as wing sections separating from each other or large kinks appearing along the span, that cause numerical deficiencies in the optimization. 
+This concern is addressed by enforcing C0 continuity along the span-wise junctions between the sections. 
+C0 continuity can be enforced for any geometric parameter OAS controls for a given section. 
+For example, the tip chord or twist of a given section should match the root chord or twist of the subsequent section. 
+There are two approaches for enforcing C0 continuity between sections: a constraint-based approach and a construction-based approach. 
+
+The constaint-based approach involves explicitly setting a position constraint that joins the surface points at section junctions to a specified tolerance. 
+This constraint is useful when varying section parameters like span or dihedral, where entire sections could collide or separate if C0 continuity is not strictly enforced. 
+A fully differentiated implementation that facilitates setting this constraint is incorporated into OAS. This approach is robust but introduces at least two linear constraints per segment junction. 
+The number of linear constraints can quickly grow for problems with many sections that this study anticipates for morphing applications. 
+In cases where geometric parameters can be specified as a function of span, however, it is possible to eliminate these additional constraints and maintain C0 continuity using the construction-based approach.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec.py
+   :start-after: checkpoint 3
+   :end-before: checkpoint 4
+
+
+The construction-based approach forgoes the constraint entirely. 
+Continuity by construction is enforced by assigning the B-spline control point located at section edges to the same independent variable controlled by the optimizer. 
+Enforcing C0 continuity by construction solely applies to geometric parameters that employ B-spline parametrization in OAS. 
+These include chord distribution, twist distribution, shear distribution in all three directions, and thickness and radius distributions for structural spars. 
+Geometric transformations parameterized by a single value, such as sweep, taper, and dihedral, cannot be used with C0 continuity enforcement by construction. 
+
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec_construction.py
+   :start-after: checkpoint 0
+   :end-before: checkpoint 1
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec_construction.py
+   :start-after: checkpoint 2
+   :end-before: checkpoint 3
+
+We can now create the aerodynamic analysis group.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec.py
+   :start-after: checkpoint 4
+   :end-before: checkpoint 5
+
+Connecting the geometry and analysis groups requires care when using the multi-section parameterization.
+While the multi-section geometry group is effectively a drop-in replacement for the standard geometery group for multi-section wings, it does require that the unified mesh component be connected to the AeroPoint analysis.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec.py
+   :start-after: checkpoint 5
+   :end-before: checkpoint 6
+
+*The following section applies to user considering cases with thickess to chord B-splines and viscous effects.*
+
+When using a thickness to chord ratio B-spline to account for viscous effects the user should be careful to connect the unified thickess to chord ratio B-spline automatically generated by the multi-section geometery group.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec_visc.py
+   :start-after: checkpoint 0
+   :end-before: checkpoint 1
+
+We can now setup our optimization problem and run it.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec.py
+   :start-after: checkpoint 6
+   :end-before: checkpoint 7
+
+We then finish by plotting the result.
+
+.. literalinclude:: /advanced_features/scripts/basic_2_sec.py
+   :start-after: checkpoint 7
+   :end-before: checkpoint 8
+
+
+
+
+
+The following shows the resulting optimized mesh. 
+The result is identical regardless of if the constraint-based or construction-based joining approahces are used.
+
+.. image:: /advanced_features/figs/multi_section_2_sym.png

openaerostruct/docs/advanced_features/scripts/basic_2_sec.py

@@ -0,0 +1,241 @@
+"""Optimizes the section chord distribution of a two section symmetrical wing using the constraint-based approach for section
+joining. This example is referenced as part of the multi-section tutorial."""
+
+# docs checkpoint 0
+import numpy as np
+import openmdao.api as om
+from openaerostruct.geometry.geometry_group import MultiSecGeometry
+from openaerostruct.aerodynamics.aero_groups import AeroPoint
+from openaerostruct.geometry.geometry_group import build_sections
+from openaerostruct.geometry.geometry_unification import unify_mesh
+import matplotlib.pyplot as plt
+
+
+# docs checkpoint 1
+
+# The multi-section geometry parameterization number section from left to right starting with section #0. A two-section symmetric wing parameterization appears as follows.
+# For a symmetrical wing the last section in the sequence will always be marked as the "root section" as it's adjacent to the geometric centerline of the wing.
+# Geometeric parameters must be specified for each section using lists with values corresponding in order of the surface numbering. Section section supports all the
+# standard OpenAeroStruct geometery transformations including B-splines.
+
+
+"""
+
+-----------------------------------------------  ^
+|                      |                       | |
+|                      |                       | |
+|        sec 0         |         sec 1         | | root         symmetrical BC
+|                      |     "root section"    | | chord
+|______________________|_______________________| |
+                                                 _
+                                              y = 0 ------------------> + y
+
+"""
+
+
+# A multi-section surface dictionary is very similar to the standard one. However, it features some additional options and requires that the user specify
+# parameters for each desired section. The multi-section geometery group also features an built in mesh generator so the wing mesh parameters can be specified right
+# in the surface dictionary. Let's create a dictionary with info and options for a two-section aerodynamic lifting surface
+
+surface = {
+    # Wing definition
+    # Basic surface parameters
+    "name": "surface",
+    "is_multi_section": True,  # This key must be present for the AeroPoint to correctly interpret this surface as multi-section
+    "num_sections": 2,  # The number of sections in the multi-section surface
+    "sec_name": [
+        "sec0",
+        "sec1",
+    ],  # names of the individual sections. Each section must be named and the list length must match the specified number of sections.
+    "symmetry": True,  # if true, model one half of wing. reflected across the midspan of the root section
+    "S_ref_type": "wetted",  # how we compute the wing area, can be 'wetted' or 'projected'
+    # Geometry Parameters
+    "taper": [1.0, 1.0],  # Wing taper for each section. The list length must match the specified number of sections.
+    "span": [2.0, 2.0],  # Wing span for each section. The list length must match the specified number of sections.
+    "sweep": [0.0, 0.0],  # Wing sweep for each section. The list length must match the specified number of sections.
+    "chord_cp": [
+        np.array([1, 1]),
+        np.array([1, 1]),
+    ],  # The chord B-spline parameterization for EACH SECTION. The list length must match the specified number of sections.
+    "twist_cp": [
+        np.zeros(2),
+        np.zeros(2),
+    ],  # The twist B-spline parameterization for EACH SECTION. The list length must match the specified number of sections.
+    "root_chord": 1.0,  # Root chord length of the section indicated as "root section"(required if using the built-in mesh generator)
+    # Mesh Parameters
+    "meshes": "gen-meshes",  # Supply a list of meshes for each section or "gen-meshes" for automatic mesh generation
+    "nx": 2,  # Number of chordwise points. Same for all sections.(required if using the built-in mesh generator)
+    "ny": [
+        21,
+        21,
+    ],  # Number of spanwise points for each section. The list length must match the specified number of sections. (required if using the built-in mesh generator)
+    # Aerodynamic Parameters
+    "CL0": 0.0,  # CL of the surface at alpha=0
+    "CD0": 0.015,  # CD of the surface at alpha=0
+    # Airfoil properties for viscous drag calculation
+    "k_lam": 0.05,  # percentage of chord with laminar
+    # flow, used for viscous drag
+    "c_max_t": 0.303,  # chordwise location of maximum (NACA0015)
+    # thickness
+    "with_viscous": False,  # if true, compute viscous drag
+    "with_wave": False,  # if true, compute wave drag
+    "groundplane": False,
+}
+
+# docs checkpoint 2
+
+# Create the OpenMDAO problem
+prob = om.Problem()
+
+# Create an independent variable component that will supply the flow
+# conditions to the problem.
+indep_var_comp = om.IndepVarComp()
+indep_var_comp.add_output("v", val=1.0, units="m/s")
+indep_var_comp.add_output("alpha", val=10.0, units="deg")
+indep_var_comp.add_output("Mach_number", val=0.3)
+indep_var_comp.add_output("re", val=1.0e5, units="1/m")
+indep_var_comp.add_output("rho", val=0.38, units="kg/m**3")
+indep_var_comp.add_output("cg", val=np.zeros((3)), units="m")
+
+# Add this IndepVarComp to the problem model
+prob.model.add_subsystem("prob_vars", indep_var_comp, promotes=["*"])
+
+# docs checkpoint 3
+
+# Instead of creating a standard geometery group, here we will create a multi-section group that will accept our multi-section surface
+# dictionary and allow us to specify any C0 continuity constraints between the sections. In this example we will constrain the sections
+# into a C0 continuous surface using a component that the optimizer can use as a constraint. The joining constraint component returns the
+# distance between the leading edge and trailing edge points at section interections. Any combination of the x,y, and z distances can be returned
+# to constrain the surface in a particular direction.
+
+
+"""
+                        LE1              LE2   cLE = [LE2x-LE1x,LE2y-LE1y,LE2z-LE1z]
+------------------------- <-------------> -------------------------
+|                       |                 |                       |
+|                       |                 |                       |
+|        sec 0          |                 |         sec 1         |
+|                       |                 |     "root section"    |
+|______________________ | <-------------> |_______________________|
+                        TE1              TE2    cTE = [TE2x-TE1x,TE2y-TE1y,TE2z-TE1z]
+
+
+
+"""
+
+# We pass in the multi-section surface dictionary to the MultiSecGeometry geometery group. We also enabled joining_comp and pass two array to dim_contr
+# These two arrays should only consists of 1 and 0 and tell the joining component which of the x,y, and z distance constraints we wish to enforce at the LE and TE
+# In this example, we only wish to constraint the x-distance between the sections at both the leading and trailing edge.
+
+multi_geom_group = MultiSecGeometry(
+    surface=surface, joining_comp=True, dim_constr=[np.array([1, 0, 0]), np.array([1, 0, 0])]
+)
+prob.model.add_subsystem(surface["name"], multi_geom_group)
+
+# docs checkpoint 4
+
+# In this next part, we will setup the aerodynamics group. First we use a utility function called build_sections which takes our multi-section surface dictionary and outputs a
+# surface dictionary for each individual section. We then inputs these dictionaries into the mesh unification function unify_mesh to produce a single mesh array for the the entire surface.
+# We then add this mesh to the multi-section surface dictionary
+section_surfaces = build_sections(surface)
+uniMesh = unify_mesh(section_surfaces)
+surface["mesh"] = uniMesh
+
+# Create the aero point group, which contains the actual aerodynamic
+# analyses. This step is exactly as it's normally done except the surface dictionary we pass in is the multi-surface one
+aero_group = AeroPoint(surfaces=[surface])
+point_name = "aero_point_0"
+prob.model.add_subsystem(point_name, aero_group, promotes_inputs=["v", "alpha", "Mach_number", "re", "rho", "cg"])
+
+# docs checkpoint 5
+
+# The following steps are similar to a normal OAS surface script but note the differences in surface naming. Note that
+# unified surface created by the multi-section geometry group needs to be connected to AeroPoint(be careful with the naming)
+
+# Get name of surface and construct the name of the unified surface mesh
+name = surface["name"]
+unification_name = "{}_unification".format(surface["name"])
+
+# Connect the mesh from the mesh unification component to the analysis point.
+prob.model.connect(name + "." + unification_name + "." + name + "_uni_mesh", point_name + "." + "surface" + ".def_mesh")
+
+# Perform the connections with the modified names within the
+# 'aero_states' group.
+prob.model.connect(
+    name + "." + unification_name + "." + name + "_uni_mesh", point_name + ".aero_states." + "surface" + "_def_mesh"
+)
+
+# docs checkpoint 6
+
+# Next, we add the DVs to the OpenMDAO problem. Note that each surface's geometeric parameters are under the given section names specified in the multi-surface dictionary earlier.
+# Here we use the chord B-spline that we specified earlier for each section and the angle-of-attack as DVs.
+prob.model.add_design_var("surface.sec0.chord_cp", lower=0.1, upper=10.0, units=None)
+prob.model.add_design_var("surface.sec1.chord_cp", lower=0.1, upper=10.0, units=None)
+prob.model.add_design_var("alpha", lower=0.0, upper=10.0, units="deg")
+
+
+# Next, we add the C0 continuity constraint for this problem by constraining the x-distance between sections to 0.
+# NOTE: SLSQP optimizer does not handle the joining equality constraint properly so the constraint needs to be specified as an inequality constraint
+# All other optimizers like SNOPT can handle the equality constraint as is.
+
+prob.model.add_constraint("surface.surface_joining.section_separation", upper=0, lower=0)  # FOR SLSQP
+# prob.model.add_constraint('surface.surface_joining.section_separation',equals=0.0,scaler=1e-4) #FOR OTHER OPTIMIZERS
+
+# Add CL constraint
+prob.model.add_constraint(point_name + ".CL", equals=0.3)
+
+# Add Wing total area constraint
+prob.model.add_constraint(point_name + ".total_perf.S_ref_total", equals=2.0)
+
+# Add objective
+prob.model.add_objective(point_name + ".CD", scaler=1e4)
+
+
+prob.driver = om.ScipyOptimizeDriver()
+prob.driver.options["optimizer"] = "SLSQP"
+prob.driver.options["tol"] = 1e-3
+prob.driver.options["disp"] = True
+prob.driver.options["maxiter"] = 1000
+prob.driver.options["debug_print"] = ["nl_cons", "objs", "desvars"]
+
+# Set up and run the optimization problem
+prob.setup()
+
+# prob.run_model()
+prob.run_driver()
+# om.n2(prob)
+
+# docs checkpoint 7
+
+# Get each section mesh
+mesh1 = prob.get_val("surface.sec0.mesh", units="m")
+mesh2 = prob.get_val("surface.sec1.mesh", units="m")
+
+# Get the unified mesh
+meshUni = prob.get_val(name + "." + unification_name + "." + name + "_uni_mesh")
+
+
+# Plot the results
+def plot_meshes(meshes):
+    """this function plots to plot the mesh"""
+    plt.figure(figsize=(8, 4))
+    for i, mesh in enumerate(meshes):
+        mesh_x = mesh[:, :, 0]
+        mesh_y = mesh[:, :, 1]
+        color = "k"
+        for i in range(mesh_x.shape[0]):
+            plt.plot(mesh_y[i, :], 1 - mesh_x[i, :], color, lw=1)
+            plt.plot(-mesh_y[i, :], 1 - mesh_x[i, :], color, lw=1)  # plots the other side of symmetric wing
+        for j in range(mesh_x.shape[1]):
+            plt.plot(mesh_y[:, j], 1 - mesh_x[:, j], color, lw=1)
+            plt.plot(-mesh_y[:, j], 1 - mesh_x[:, j], color, lw=1)  # plots the other side of symmetric wing
+    plt.axis("equal")
+    plt.xlabel("y (m)")
+    plt.ylabel("x (m)")
+    plt.savefig("opt_planform_constraint.png")
+
+
+# plot_meshes([mesh1,mesh2])
+plot_meshes([meshUni])
+# plt.show()
+# docs checkpoint 8