Source code for openaerostruct.structures.wing_weight_loads

import numpy as np

from scipy.sparse import coo_matrix, diags

import openmdao.api as om
from openaerostruct.structures.utils import norm
from openaerostruct.utils.constants import grav_constant




class StructureWeightLoads(om.ExplicitComponent):
    """
    Compute the nodal loads from the weight of the wing structure to be applied to the wing
    structure.

    Parameters
    ----------
    element_mass[ny-1] : numpy array
        Weight for each wing-structure segment.
    nodes[ny, 3] : numpy array
        Flattened array with coordinates for each FEM node.
    load_factor : float
        Load factor for the flight point.

    Returns
    -------
    struct_weight_loads[ny, 6] : numpy array
        Flattened array containing the loads applied on the FEM component,
        computed from the weight of the wing-structure segments.
    element_lengths[ny-1] : numpy array
        Lengths of the FEM finite elements.
    """

    def initialize(self):
        self.options.declare("surface", types=dict)

    def setup(self):
        self.surface = surface = self.options["surface"]
        self.ny = surface["mesh"].shape[1]

        self.add_input("element_mass", val=np.zeros((self.ny - 1)), units="kg")
        self.add_input("nodes", val=np.zeros((self.ny, 3), dtype=complex), units="m")
        self.add_input("load_factor", val=1.0)

        self.add_output("struct_weight_loads", val=np.zeros((self.ny, 6)), units="N")
        self.add_output("element_lengths", val=np.zeros(self.ny - 1), units="m")

        nym1 = self.ny - 1

        # dloads__dzf
        rows = np.zeros(2 * nym1)
        cols = np.zeros(2 * nym1)
        data = -np.ones(2 * nym1)
        counter = 0
        for i in range(nym1):
            for j in range(2):
                rows[counter] = 6 * (i + j) + 2
                cols[counter] = i
                counter += 1

        self.dswl__dzf = coo_matrix((data, (rows, cols)), shape=(self.ny * 6, nym1))

        rows = np.zeros(3 * self.ny)
        cols = np.zeros(3 * self.ny)
        for j in range(3):
            for i in range(self.ny):
                idx = i + self.ny * j
                rows[idx] = 6 * i + j + 2

        self.dswl__dlf_row = rows
        self.dswl__dlf_col = cols

        self.declare_partials("struct_weight_loads", "load_factor", rows=rows, cols=cols)
        # self.declare_partials('struct_weight_loads', 'load_factor')

        # dstruct_weight_loads__dnodes (this one is super complicated)
        #     we can pre-compute several of the intermediate partial derivative matrices
        #     others we can pre-compute the rows/cols so we just have to set the data
        #     the for loops for rows/cols are a bit expensive (cause python), but the data computations are fast
        #     so we will pay the one-time setup cost to pre-compute the rows/cols
        rows_del = np.zeros(2 * nym1)
        cols_d0 = np.zeros(2 * nym1)
        cols_d1 = np.zeros(2 * nym1)
        data = np.zeros(2 * nym1)

        counter = 0
        for i in range(nym1):
            for j in range(2):
                rows_del[counter] = i
                cols_d0[counter] = 3 * (i + j)
                cols_d1[counter] = 3 * (i + j) + 1
                data[counter] = 2 * j - 1  # this gives 1 or -1 alternating based on j
                counter += 1

        shape = (nym1, 3 * self.ny)
        self.ddel0__dnodes = coo_matrix((data, (rows_del, cols_d0)), shape=shape)
        self.ddel1__dnodes = coo_matrix((data, (rows_del, cols_d1)), shape=shape)

        # del__dnodes
        rows_el = np.zeros(6 * nym1)
        cols_el = np.zeros(6 * nym1)
        data = np.ones(6 * nym1)  # del__dnodes matrix has only ones in it

        counter = 0
        for i in range(nym1):
            for j in range(6):
                rows_el[counter] = i
                cols_el[counter] = 3 * i + j
                counter += 1
        self.del__dnodes = coo_matrix((data, (rows_el, cols_el)), shape=shape)

        rng = np.random.default_rng(0)
        nym1_rand = rng.random(nym1)

        dzm_dnodes_pattern = (
            (diags(nym1_rand * nym1_rand) * self.ddel0__dnodes + diags(nym1_rand * nym1_rand) * self.ddel1__dnodes)
        ).tocsr()

        dbm3_dnodes = (
            diags(nym1_rand) * dzm_dnodes_pattern
            + diags(nym1_rand) * self.ddel1__dnodes
            - diags(nym1_rand) * self.del__dnodes
        )

        dbm4_dnodes = (
            diags(nym1_rand) * dzm_dnodes_pattern
            + diags(nym1_rand) * self.ddel0__dnodes
            - diags(nym1_rand) * self.del__dnodes
        )

        self.dbm3_dnodes_pattern = (
            diags(nym1_rand) * dzm_dnodes_pattern
            + diags(nym1_rand) * self.ddel1__dnodes
            - diags(nym1_rand) * self.del__dnodes
        )
        self.dbm4_dnodes_pattern = (
            diags(nym1_rand) * dzm_dnodes_pattern
            + diags(nym1_rand) * self.ddel0__dnodes
            - diags(nym1_rand) * self.del__dnodes
        )
        # dswl__dbm
        rows_1 = np.zeros(2 * nym1)
        rows_0 = np.zeros(2 * nym1)
        cols = np.zeros(2 * nym1)
        data = rng.random(2 * nym1)
        counter = 0
        for i in range(nym1):
            for j in range(2):
                rows_1[counter] = 6 * (i + j) + 3
                rows_0[counter] = 6 * (i + j) + 4
                if counter % 2 == 0:  # even
                    data[counter] = 1
                else:  # odd
                    data[counter] = -1
                cols[counter] = i
                counter += 1

        self.dswl__dbm3 = coo_matrix((data, (rows_1, cols)), shape=(6 * self.ny, nym1)).tocsr()
        self.dswl__dbm4 = coo_matrix((data, (rows_0, cols)), shape=(6 * self.ny, nym1)).tocsr()

        self.dswl__dnodes_pattern = (self.dswl__dbm3 * dbm3_dnodes + self.dswl__dbm4 * dbm4_dnodes).tocoo()

        self.dswl__dew_pattern = (self.dswl__dbm4 + self.dswl__dbm3 + self.dswl__dzf).tocoo()

        self.declare_partials(
            "struct_weight_loads", "nodes", rows=self.dswl__dnodes_pattern.row, cols=self.dswl__dnodes_pattern.col
        )

        self.declare_partials(
            "struct_weight_loads", "element_mass", rows=self.dswl__dew_pattern.row, cols=self.dswl__dew_pattern.col
        )
        self.set_check_partial_options(wrt="*", method="cs")

    def compute(self, inputs, outputs):
        struct_weights = inputs["element_mass"] * inputs["load_factor"] * grav_constant
        nodes = inputs["nodes"]

        element_lengths = norm(nodes[1:, :] - nodes[:-1, :], axis=1)

        # And we also need the deltas between consecutive nodes
        deltas = nodes[1:, :] - nodes[:-1, :]
        # save these slices cause I use them a lot
        del0 = deltas[:, 0]
        del1 = deltas[:, 1]

        # Assume weight coincides with the elastic axis
        z_forces_for_each = struct_weights / 2.0
        z_moments_for_each = struct_weights / 12.0 * (del0**2 + del1**2) ** 0.5

        loads = outputs["struct_weight_loads"]
        loads *= 0  # need to zero it out, since we're accumulating onto it
        # Why doesn't this trigger when running ../../tests/test_aerostruct_wingbox_+weight_analysis.py???
        if self.under_complex_step:
            loads = np.zeros((self.ny, 6), dtype=complex)

        # Loads in z-direction
        loads[:-1, 2] += -z_forces_for_each
        loads[1:, 2] += -z_forces_for_each

        # Bending moments for consistency
        bm3 = z_moments_for_each * del1 / element_lengths
        loads[:-1, 3] += -bm3
        loads[1:, 3] += bm3

        bm4 = z_moments_for_each * del0 / element_lengths
        loads[:-1, 4] += -bm4
        loads[1:, 4] += bm4
        outputs["struct_weight_loads"] = loads

    def compute_partials(self, inputs, J):
        struct_weights = inputs["element_mass"] * inputs["load_factor"] * grav_constant
        nodes = inputs["nodes"]

        element_lengths = norm(nodes[1:, :] - nodes[:-1, :], axis=1)

        # And we also need the deltas between consecutive nodes
        deltas = nodes[1:, :] - nodes[:-1, :]
        # save these slices cause I use them alot
        del0 = deltas[:, 0]
        del1 = deltas[:, 1]

        # Assume weight coincides with the elastic axis
        z_moments_for_each = struct_weights / 12.0 * (del0**2 + del1**2) ** 0.5

        dzf__dew = 0.5 * inputs["load_factor"][0]
        dzf__dlf = inputs["element_mass"] / 2.0

        dzm__dew = diags((del0**2 + del1**2) ** 0.5 / 12 * inputs["load_factor"])
        dzm__dlf = (del0**2 + del1**2) ** 0.5 / 12.0 * inputs["element_mass"]

        dbm3__dzm = diags(del1 / element_lengths)
        dbm4__dzm = diags(del0 / element_lengths)

        # need to convert to lil to re-order the data to match the original row/col indexing from setup
        dswl__dlf = (
            -self.dswl__dbm3 * dbm3__dzm * coo_matrix(dzm__dlf).T
            + -self.dswl__dbm4 * dbm4__dzm * coo_matrix(dzm__dlf).T
            + self.dswl__dzf * coo_matrix(dzf__dlf).T
        ).tolil()

        J["struct_weight_loads", "load_factor"] = (
            dswl__dlf[self.dswl__dlf_row, self.dswl__dlf_col].toarray().flatten() * grav_constant
        )

        dswl__dew = (
            -self.dswl__dbm4 * dbm4__dzm * dzm__dew
            + -self.dswl__dbm3 * dbm3__dzm * dzm__dew
            + self.dswl__dzf * dzf__dew
        ).tolil()

        data = dswl__dew[self.dswl__dew_pattern.row, self.dswl__dew_pattern.col].toarray().flatten()
        J["struct_weight_loads", "element_mass"] = data * grav_constant

        # dstruct_weight_loads__dnodes (this one is super complicated)

        # del__dnodes
        # note: del__dnodes matrix already created in setup, just need to set data
        raw_data = diags(1 / element_lengths) * (nodes[1:, :] - nodes[:-1, :])
        data = np.hstack((-raw_data, raw_data)).flatten()
        self.del__dnodes.data = data

        dzm_dnodes = diags(struct_weights / 12 * (del0**2 + del1**2) ** -0.5) * (
            diags(del0) * self.ddel0__dnodes + diags(del1) * self.ddel1__dnodes
        )

        dbm3_dnodes = (
            diags(del1 / element_lengths) * dzm_dnodes
            + diags(z_moments_for_each / element_lengths) * self.ddel1__dnodes
            - diags(z_moments_for_each * del1 / element_lengths**2) * self.del__dnodes
        )

        dbm4_dnodes = (
            diags(del0 / element_lengths) * dzm_dnodes
            + diags(z_moments_for_each / element_lengths) * self.ddel0__dnodes
            - diags(z_moments_for_each * del0 / element_lengths**2) * self.del__dnodes
        )

        # this is kind of dumb, but I need lil cause I have to re-index to preserve order
        #     the coo column ordering doesn't seem to be deterministic
        #     so I use the original row/col from the pattern as index arrays to
        #     pull the data out in the correct order
        dswl__dnodes = (self.dswl__dbm3 * dbm3_dnodes + self.dswl__dbm4 * dbm4_dnodes).tolil()
        data = dswl__dnodes[self.dswl__dnodes_pattern.row, self.dswl__dnodes_pattern.col].toarray().flatten()
        J["struct_weight_loads", "nodes"] = -data