From a4a77d9eab760ee9ca717e688e6bd43990ed1ec6 Mon Sep 17 00:00:00 2001
From: mjburton <mjburton@mit.edu>
Date: Thu, 2 Nov 2017 19:47:48 -0400
Subject: [PATCH] modularize beam model; to allow for different qbar and Sbar
 functions

---
 .../GP/aircraft/wing/chord_spar_loading.py    | 17 ++++++++++++--
 gpkitmodels/GP/aircraft/wing/wing_test.py     |  8 ++++---
 gpkitmodels/GP/beam/beam.py                   | 22 ++++++++++++++-----
 3 files changed, 36 insertions(+), 11 deletions(-)

diff --git a/gpkitmodels/GP/aircraft/wing/chord_spar_loading.py b/gpkitmodels/GP/aircraft/wing/chord_spar_loading.py
index d89cafb2..ba65d157 100644
--- a/gpkitmodels/GP/aircraft/wing/chord_spar_loading.py
+++ b/gpkitmodels/GP/aircraft/wing/chord_spar_loading.py
@@ -1,4 +1,5 @@
 " cap spar "
+from numpy import flip
 from gpkit import Model, Variable, Vectorize
 from gpkitmodels.GP.beam.beam import Beam
 
@@ -17,12 +18,24 @@ def setup(self, static, Wcent):
         with Vectorize(static.N-1):
             Mr = Variable("M_r", "N*m", "wing section root moment")
 
+        def new_qbarFun(_, c):
+            " define qbar model for chord loading "
+            return [f(c) for f in static.substitutions["\\bar{c}"]]
+        def new_SbarFun(bmodel, c):
+            " define Sbar model for chord loading "
+            Sb = [1e-10]*static.N
+            for i in flip(range(static.N-1), 0):
+                Sb[i] = (Sb[i+1] + static.substitutions["d\\eta"][i](c)
+                         * (new_qbarFun(bmodel, c)[i]
+                            + new_qbarFun(bmodel, c)[i+1])/2)
+            return Sb
+
+        Beam.qbarFun = new_qbarFun
+        Beam.SbarFun = new_SbarFun
         beam = Beam(static.N)
-        beam.substitutions["\\bar{q}"] = static.substitutions["\\bar{c}"]
 
         constraints = [
             # dimensionalize moment of inertia and young's modulus
-            # beam["\\bar{q}"] == static["\\bar{c}"],
             beam["dx"] == static["d\\eta"],
             beam["\\bar{EI}"] <= (8*static["E"]*static["I"]/Nmax
                                   / Wcent/static["b"]**2),
diff --git a/gpkitmodels/GP/aircraft/wing/wing_test.py b/gpkitmodels/GP/aircraft/wing/wing_test.py
index 3fbe9576..b5339524 100644
--- a/gpkitmodels/GP/aircraft/wing/wing_test.py
+++ b/gpkitmodels/GP/aircraft/wing/wing_test.py
@@ -23,10 +23,12 @@ def wing_test():
     W.substitutions[W.topvar("W")] = 50
     fs = FlightState()
     perf = W.flight_model(W, fs)
-    loading = W.loading(W, Wcent, W.topvar("W"), fs["V"], perf["C_L"])
+    # loading = W.loading(W, Wcent, W.topvar("W"), fs["V"], perf["C_L"])
+    loading = W.loading(W, Wcent)
 
     m = Model(perf["C_d"], [W, fs, perf, loading])
-    m.solve("mosek")
+    s = m.solve("mosek")
+    return s
 
 if __name__ == "__main__":
-    wing_test()
+    s = wing_test()
diff --git a/gpkitmodels/GP/beam/beam.py b/gpkitmodels/GP/beam/beam.py
index 40ad0c94..34bd2546 100644
--- a/gpkitmodels/GP/beam/beam.py
+++ b/gpkitmodels/GP/beam/beam.py
@@ -1,8 +1,13 @@
 " discretized beam model "
 from gpkit import Model, Variable, Vectorize
 
+#pylint: disable=invalid-name
+
 class Beam(Model):
     "discretized beam bending model"
+    qbarFun = None
+    SbarFun = None
+
     def setup(self, N):
 
         with Vectorize(N-1):
@@ -11,8 +16,8 @@ def setup(self, N):
             dx = Variable("dx", "-", "normalized length of element")
 
         with Vectorize(N):
-            qbar = Variable("\\bar{q}", "-", "normalized loading")
-            Sbar = Variable("\\bar{S}", "-", "normalized shear")
+            qbar = Variable("\\bar{q}", self.qbarFun, "-", "normalized loading")
+            Sbar = Variable("\\bar{S}", self.SbarFun, "-", "normalized shear")
             Mbar = Variable("\\bar{M}", "-", "normalized moment")
             th = Variable("\\theta", "-", "deflection slope")
             dbar = Variable("\\bar{\\delta}", "-", "normalized displacement")
@@ -24,15 +29,20 @@ def setup(self, N):
         dbarroot = Variable("\\bar{\\delta}_{root}", 1e-10, "-",
                             "Base deflection")
 
-        constraints = [
-            Sbar[:-1] >= Sbar[1:] + 0.5*dx*(qbar[:-1] + qbar[1:]),
-            Sbar[-1] >= Sbartip,
+        constraints = []
+        if self.SbarFun is None:
+            constraints.extend([
+                Sbar[:-1] >= Sbar[1:] + 0.5*dx*(qbar[:-1] + qbar[1:]),
+                Sbar[-1] >= Sbartip])
+
+        constraints.extend([
             Mbar[:-1] >= Mbar[1:] + 0.5*dx*(Sbar[:-1] + Sbar[1:]),
             Mbar[-1] >= Mbartip,
             th[0] >= throot,
             th[1:] >= th[:-1] + 0.5*dx*(Mbar[1:] + Mbar[:-1])/EIbar,
             dbar[0] >= dbarroot,
             dbar[1:] >= dbar[:-1] + 0.5*dx*(th[1:] + th[:-1]),
-            ]
+            ])
 
         return constraints
+