Merge branch 'main' into joss-paper

.github/test_real.sh

@@ -2,8 +2,8 @@
 set -e
 ./input_files/get-input-files.sh
 
-# No tests should be skipped on GCC and non Intel MPI
-if [[ $COMPILERS == "gcc" ]] && [[ -z $I_MPI_ROOT ]]; then
+# No tests should be skipped on GCC, non Intel MPI, and x86 arch
+if [[ $COMPILERS == "gcc" ]] && [[ -z $I_MPI_ROOT ]] && [[ "$(arch)" == "x86_64" ]]; then
     EXTRA_FLAGS='--disallow_skipped'
 fi
 

.gitignore

@@ -6,3 +6,4 @@ reg_tests/pygeo_reg.orig
 doc/_build
 *testflo_report.out
 input_files
+.isort.cfg

doc/advanced_ffd.rst

@@ -15,7 +15,7 @@ During this tutorial, we will use the Cessna 172 airplane wing as our example ge
 .. figure:: images/c172.jpg
    :width: 450
    :target: images/c172.jpg
-   :align: center 
+   :align: center
 
    A Cessna 172 in flight (Anna Zvereva, CC-BY)
 
@@ -27,7 +27,7 @@ We use the spline surfacing tool in `pyGeo` to generate the Cessna wing geometry
 A full description of the surfacing script is beyond the scope of this tutorial, but the script itself can be found at ``examples/c172_wing/c172.py``.
 Using the ``pyiges`` package, one of the IGES CAD files can be converted to a triangulated (.stl) format in order to turn the wing into a pointset.
 
-Next, we need to create an FFD volume that encloses the wing. 
+Next, we need to create an FFD volume that encloses the wing.
 We want to approximate the wing closely without any of the wing intersecting the box.
 Using our knowledge of the wing dimensions, it's easy to create a closely-conforming FFD.
 The example script is located at ``examples/c172_wing/genFFD.py``.
@@ -56,7 +56,7 @@ In the following example, we perturb a single point in the inner portion of the
 
 .. literalinclude:: ../examples/c172_wing/runFFDExample.py
     :start-after: # rst add local DV
-    :end-before: # rst ref axis
+    :end-before: # rst add shape function DVs
 
 
 This local perturbation produces the obvious deformation in the following rendering:
@@ -68,6 +68,17 @@ This local perturbation produces the obvious deformation in the following render
 
 .. _global_vars:
 
+--------------------------------------
+Adding shape function design variables
+--------------------------------------
+
+Similar to local design variables, shape function DVs can be used to define design variables that control the local shape using a single or several control points in specified directions. The shape functions define a direction and magnitude for control points. This DV can be used to link control point movements without using linear constraints.
+
+.. literalinclude:: ../examples/c172_wing/runFFDExample.py
+    :start-after: # rst add shape function DVs
+    :end-before: # rst ref axis
+
+
 -----------------------------------
 Reference axes and global variables
 -----------------------------------
@@ -83,7 +94,7 @@ Global design variables commonly include the following mathematical transformati
 - Linear stretching or shrinking
 - Translation
 
-Rotating a point requires knowing an axis of rotation. 
+Rotating a point requires knowing an axis of rotation.
 Scaling a point requires a reference point.
 We can define these for the entire pointset by defining one or more *reference axes*.
 A reference axis is defined as a line or curve within the FFD volume.
@@ -142,7 +153,7 @@ The global design variable can be perturbed just like a local design variable, a
     :start-after: # rst set DV
     :end-before: # rst set DV 2
 
-Applying this twist results in the geometry pictured below. 
+Applying this twist results in the geometry pictured below.
 The location of the reference axis (and any points located close to the reference axis) is not affected by the rotation.
 This is a general principle of applying transformations: *the reference axis location remains invariant under the transformation*.
 
@@ -192,12 +203,12 @@ The results of the sweep are dramatic, as seen in the rendering.
 This example illustrates an important detail; namely, that the local control points do not rotate in the x-z plane as the wing is swept back.
 This is because of the way the reference axis is implemented.
 Every local control point (the red dots) is *projected* onto the reference axis when the axis is created.
-In this case, by default, the points were projected along the x axis. 
+In this case, by default, the points were projected along the x axis.
 Once the points are projected, they become rigidly linked to the projected point on the axis.
 Even if the reference axis is rotated, the rigid links do *not* rotate.
 However, the links do translate along with their reference point.
 Only the ``scale_`` and ``rot_`` operators change the rigid links.
-:meth:`writeLinks <.DVGeometry.writeLinks>` can be used to write out these links, which can then be viewed in Tecplot. 
+:meth:`writeLinks <.DVGeometry.writeLinks>` can be used to write out these links, which can then be viewed in Tecplot.
 
 --------------------------------
 Multiple global design variables
@@ -221,7 +232,7 @@ Let's also introduce a random perturbation to the local design variables to see
 .. literalinclude:: ../examples/c172_wing/runFFDExample.py
     :start-after: # rst set DV 4
 
-The combination of multiple global and local design variables produces the wild shape in the rendering below. 
+The combination of multiple global and local design variables produces the wild shape in the rendering below.
 Obviously this is not a suitable optimized aircraft design.
 However, the optimizer is free to use all of these degrees of freedom to eventually find the best possible result.
 

examples/c172_wing/runFFDExample.py

@@ -46,6 +46,75 @@ def create_fresh_dvgeo():
 stlmesh.save("local_wing.stl")
 DVGeo.writeTecplot("local_ffd.dat")
 
+
+# rst add shape function DVs
+DVGeo, stlmesh = create_fresh_dvgeo()
+# create the shape functions. to demonstrate the capability, we add 2 shape functions. each
+# adds a bump to a spanwise section, and then the two neighboring spanwise sections also get
+# half of the perturbation.
+
+# this array can be access as [i,j,k] that follows the FFD volume's topology, and returns the global
+# coef indices, which is what we need to set the shape
+lidx = DVGeo.getLocalIndex(0)
+
+# create the dictionaries
+shape_1 = {}
+shape_2 = {}
+
+k_center = 4
+i_center = 5
+n_chord = lidx.shape[0]
+
+for kk in [-1, 0, 1]:
+    if kk == 0:
+        k_weight = 1.0
+    else:
+        k_weight = 0.5
+
+    for ii in range(n_chord):
+        # compute the chord weight. we want the shape to peak at i_center
+        if ii == i_center:
+            i_weight = 1.0
+        elif ii < i_center:
+            # we are ahead of the center point
+            i_weight = ii / i_center
+        else:
+            # we are behind the center point
+            i_weight = (n_chord - ii - 1) / (n_chord - i_center - 1)
+
+        # get the direction vectors with unit length
+        dir_up = np.array([0.0, 1.0, 0.0])
+        # dir down can also be defined as an upwards pointing vector. Then, the DV itself
+        # getting a negative value means the surface would move down etc. For now, we define
+        # the vector as its pointing down, so a positive DV value moves the surface down.
+        dir_down = np.array([0.0, -1.0, 0.0])
+
+        # scale them by the i and k weights
+        dir_up *= k_weight * i_weight
+        dir_down *= k_weight * i_weight
+
+        # get this point's global index and add to the dictionary with the direction vector.
+        gidx_up = lidx[ii, 1, kk + k_center]
+        gidx_down = lidx[ii, 0, kk + k_center]
+
+        shape_1[gidx_up] = dir_up
+        # the lower face is perturbed with a separate dictionary
+        shape_2[gidx_down] = dir_down
+
+shapes = [shape_1, shape_2]
+DVGeo.addShapeFunctionDV("shape_func", shapes)
+
+dvdict = DVGeo.getValues()
+dvdict["shape_func"] = np.array([0.3, 0.2])
+DVGeo.setDesignVars(dvdict)
+
+# write out to data files for visualization
+stlmesh.vectors[:, 0, :] = DVGeo.update("mesh_v0")
+stlmesh.vectors[:, 1, :] = DVGeo.update("mesh_v1")
+stlmesh.vectors[:, 2, :] = DVGeo.update("mesh_v2")
+stlmesh.save("shape_func_wing.stl")
+DVGeo.writeTecplot("shape_func_ffd.dat")
+
 # rst ref axis
 DVGeo, stlmesh = create_fresh_dvgeo()
 # add a reference axis named 'c4' to the FFD volume

pygeo/mphys/mphys_dvgeo.py

@@ -310,6 +310,50 @@ def nom_addLocalSectionDV(
         self.add_input(dvName, distributed=False, shape=nVal)
         return nVal
 
+    def nom_addShapeFunctionDV(self, dvName, shapes, childIdx=None, config=None):
+        """
+        Add one or more local shape function design variables to the DVGeometry object
+        Wrapper for :meth:`addShapeFunctionDV <.DVGeometry.addShapeFunctionDV>`
+        Input parameters are identical to those in wrapped function unless otherwise specified
+
+        Parameters
+        ----------
+        dvName : str
+            Name to give this design variable
+        shapes : list of dictionaries, or a single dictionary
+            See wrapped
+        childIdx : int, optional
+            The zero-based index of the child FFD, if this DV is for a child FFD
+            The index is defined by the order in which you add the child FFD to the parent
+            For example, the first child FFD has an index of 0, the second an index of 1, and so on
+        config : str or list, optional
+            See wrapped
+
+        Returns
+        -------
+        N : int
+            The number of design variables added.
+
+        Raises
+        ------
+        RuntimeError
+            Raised if the underlying DVGeo parameterization is not FFD-based
+        """
+        # shape function DVs are only added to FFD-based DVGeo objects
+        if self.geo_type != "ffd":
+            raise RuntimeError(f"Only FFD-based DVGeo objects can use local DVs, not type:{self.geo_type}")
+
+        # add the DV to a normal DVGeo
+        if childIdx is None:
+            nVal = self.DVGeo.addShapeFunctionDV(dvName, shapes, config)
+        # add the DV to a child DVGeo
+        else:
+            nVal = self.DVGeo.children[childIdx].addShapeFunctionDV(dvName, shapes, config)
+
+        # define the input
+        self.add_input(dvName, distributed=False, shape=nVal)
+        return nVal
+
     def nom_addGeoCompositeDV(self, dvName, ptSetName=None, u=None, scale=None, **kwargs):
         # call the dvgeo object and add this dv
         self.DVGeo.addCompositeDV(dvName, ptSetName=ptSetName, u=u, scale=scale, **kwargs)