Merge pull request #8 from farg-bristol/develop

pull in v0.1.5

.github/workflows/python-publish.yml

@@ -0,0 +1,32 @@
+# This workflow will upload a Python Package using Twine when a release is created
+# For more information see: https://help.github.com/en/actions/language-and-framework-guides/using-python-with-github-actions#publishing-to-package-registries
+
+name: Upload Python Package
+
+on:
+  release:
+    types: [created]
+
+jobs:
+  deploy:
+
+    runs-on: ubuntu-latest
+
+    steps:
+    - uses: actions/checkout@v2
+    - name: Set up Python
+      uses: actions/setup-python@v2
+      with:
+        python-version: '3.x'
+    - name: Install dependencies
+      run: |
+        python -m pip install --upgrade pip
+        pip install setuptools wheel twine
+    - name: Build package
+      run: |
+        python setup.py sdist bdist_wheel
+    - name: Publish package
+      uses: pypa/gh-action-pypi-publish@release/v1
+      with:
+        user: __token__
+        password: ${{ secrets.PYPI_API_TOKEN }}

changelog.txt

@@ -0,0 +1,32 @@
+# Changelog
+
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to [Explicit Versioning](https://github.com/exadra37-versioning/explicit-versioning).
+
+## [0.1.5.0]
+### Added
+- added optional flag on the 'BodyForce' constructor to coose whether to 
+    simplify the resulting expresion
+### Changed
+- updated octave / python print methods to utilise common sub expersion 
+    decompositon. This significantly increase the performance of the generated 
+    EoMs
+### Deprecated
+### Removed
+### Fixed
+
+## [0.1.4.2] - 2021-04-07
+### Added
+- first commit to moyra: this software was created to easily create the EoMs of
+    basic nonlinear aeroelastic systems using the symbolic library sympy
+- added methods the create multibody elements and forces and create the EoMs via
+    the Euler-Lagrange method
+- support to generate forces from basic strip theory
+- ability to print the generated symbolic EoMs to either numpy or MATLAB/Octave 
+    to use the EoM with numerical methods
+### Changed
+### Deprecated
+### Removed
+### Fixed

moyra/__init__.py

@@ -1,4 +1,4 @@
-from .model_parameters import ModelParameters, ModelMatrix, ModelSymbol,ModelMatrixSymbol, ModelExpr
+from .model_parameters import ModelParameters, ModelMatrix, ModelSymbol,ModelMatrixSymbol, ModelExpr, ModelValue
 from .dynamic_model_parameters import DynamicModelParameters
 from .helper_funcs import linearise_matrix, extract_eigen_value_data
 from .symbolic_model import SymbolicModel

moyra/dynamic_model_parameters.py

@@ -34,10 +34,10 @@ def print_python(self):
                     else:
                         string += f'\np.{name} = ma.ModelMatrix(value={value.value}, string=\'{value._matrix_symbol}\', length={len(value)})'
                 elif isinstance(value,sym.Symbol):
-                    string += f'\np.{name} = sym.Symbol(\'{value.name}\')'
+                    string += f'\np.{name} = Symbol(\'{value.name}\')'
         return string
 
     def to_file(self,filename,file_dir=''):
-        string = 'import moyra as ma\n\ndef get_p():\n\t' + self.print_python().replace('\n','\n\t') + '\n\treturn p\n'
+        string = 'import moyra as ma\nfrom sympy import *\n\ndef get_p():\n\t' + self.print_python().replace('\n','\n\t') + '\n\treturn p\n'
         with open(os.path.join(file_dir,filename),'w') as file:
             file.write(string)

moyra/forces/aero_force.py

@@ -34,12 +34,46 @@ def PerUnitSpan(cls,p,Transform,C_L,alphadot,M_thetadot,e,rootAlpha,alpha_zero =
         """
         ## force per unit length will following theredosons pseado-steady theory
 
+        BodyJacobian = sym.simplify(cls._trigsimp(Transform.BodyJacobian(p.q)))
+
+        (wrench,dAlpha) = cls.get_wrench(p,BodyJacobian,C_L,alphadot,M_thetadot,e,
+                            rootAlpha,alpha_zero,stall_angle,c_d_max,w_g,
+                            V,c,linear,z_inverted)
+
+        _Q = BodyJacobian.T*wrench
+
+        return cls(_Q,dAlpha)
+
+    def __init__(self,Q,dAlpha):
+        self.dAlpha = dAlpha
+        super().__init__(Q) 
+
+    @staticmethod
+    def _trigsimp(expr):
+        return sym.trigsimp(sym.powsimp(sym.cancel(sym.expand(expr))))      
+
+    def linearise(self,x,x_f):
+        Q_lin = linearise_matrix(self.Q(),x,x_f)
+        dAlpha_lin = linearise_matrix(self.dAlpha,x,x_f)
+        return AeroForce(Q_lin,dAlpha_lin)
+
+    def subs(self,*args):
+        return AeroForce(self._Q.subs(*args),self.dAlpha.subs(*args))
+
+    def msubs(self,*args):
+        return AeroForce(me.msubs(self._Q,*args),me.msubs(self.dAlpha,*args))
+
+    def integrate(self,*args):
+        return AeroForce(self._Q.integrate(*args),self.dAlpha)
+
+    @staticmethod
+    def get_wrench(p,BodyJacobian,C_L,alphadot,M_thetadot,e,rootAlpha,alpha_zero = 0,stall_angle=0.24,c_d_max = 1,w_g = 0,V=None,c=None,linear=False,z_inverted=False):
+        """
+        see the class method PerUnitSpan for a explaination of terms
+        """
         if c is None:
             c=p.c
 
-        # add z velocity due to motion
-        BodyJacobian = sym.simplify(cls._trigsimp(Transform.BodyJacobian(p.q)))
-
         v_z_eff = (BodyJacobian*p.qd)[2]
         if z_inverted:
             v_z_eff *= -1
@@ -58,13 +92,11 @@ def PerUnitSpan(cls,p,Transform,C_L,alphadot,M_thetadot,e,rootAlpha,alpha_zero =
         else:
             c_l = C_L*(1/p.clip_factor*sym.ln((1+sym.exp(p.clip_factor*(dAlpha+stall_angle)))/(1+sym.exp(p.clip_factor*(dAlpha-stall_angle))))-stall_angle)
 
-        c_d = c_d_max*sym.Rational(1,2)*(1-sym.cos(2*dAlpha))
-
-        ang = rootAlpha - v_z_eff/V
-
         if linear:
             c_n = c_l
         else:
+            c_d = c_d_max*sym.Rational(1,2)*(1-sym.cos(2*dAlpha))
+            ang = rootAlpha - v_z_eff/V
             c_n = c_l*sym.cos(ang)+c_d*sym.sin(ang)
 
         F_n = dynamicPressure*c*c_n
@@ -77,32 +109,8 @@ def PerUnitSpan(cls,p,Transform,C_L,alphadot,M_thetadot,e,rootAlpha,alpha_zero =
             wrench = sym.Matrix([0,0,-F_n,0,M_w,0])
         else:
             wrench = sym.Matrix([0,0,F_n,0,M_w,0])
+        return (wrench,dAlpha)
 
-        _Q = BodyJacobian.T*wrench
-
-        return cls(_Q,dAlpha)
-
-    def __init__(self,Q,dAlpha):
-        self.dAlpha = dAlpha
-        super().__init__(Q) 
-
-    @staticmethod
-    def _trigsimp(expr):
-        return sym.trigsimp(sym.powsimp(sym.cancel(sym.expand(expr))))      
-
-    def linearise(self,x,x_f):
-        Q_lin = linearise_matrix(self.Q(),x,x_f)
-        dAlpha_lin = linearise_matrix(self.dAlpha,x,x_f)
-        return AeroForce(Q_lin,dAlpha_lin)
-
-    def subs(self,*args):
-        return AeroForce(self._Q.subs(*args),self.dAlpha.subs(*args))
-
-    def msubs(self,*args):
-        return AeroForce(me.msubs(self._Q,*args),me.msubs(self.dAlpha,*args))
-
-    def integrate(self,*args):
-        return AeroForce(self._Q.integrate(*args),self.dAlpha)
 
 
 

moyra/forces/body_force.py

@@ -5,7 +5,7 @@
 
 class BodyForce(ExternalForce):
     """A class used to represent Forces and moment in a particular reference frame"""        
-    def __init__(self,p,Transform,Fx=0,Fy=0,Fz=0,Mx=0,My=0,Mz=0):
+    def __init__(self,p,Transform,Fx=0,Fy=0,Fz=0,Mx=0,My=0,Mz=0,simplify=True):
         """
         Constructor for a body force, with the following parameters:
 
@@ -20,7 +20,11 @@ def __init__(self,p,Transform,Fx=0,Fy=0,Fz=0,Mx=0,My=0,Mz=0):
         My          - (default = 0) the moment applied to the body about the local y axis
         Mz          - (default = 0) the moment applied to the body about the local z axis
         """
-        BodyJacobian = sym.simplify(self._trigsimp(Transform.BodyJacobian(p.q)))
+
+        if simplify:
+            BodyJacobian = sym.simplify(self._trigsimp(Transform.BodyJacobian(p.q)))
+        else:
+            BodyJacobian = Transform.BodyJacobian(p.q)
         wrench = sym.Matrix([Fx,Fy,Fz,Mx,My,Mz])
         super().__init__(BodyJacobian.T*wrench) 
 

moyra/model_parameters.py

@@ -36,7 +36,7 @@ def __eq__(self,other):
     def __hash__(self):
         return hash(sym.Symbol(self.name))
     def _octave(self,printer):
-        return f'p.{self.name}'
+        return f'{self.name}'
 
 class ModelMatrixSymbol(ModelSymbol):
     def __init__(self,string,**kwarg):
@@ -46,7 +46,7 @@ def __init__(self,string,**kwarg):
     def __new__(cls,string,**kwarg):
         return super().__new__(cls,string,**kwarg)
     def _octave(self,printer):
-        return f'p.{self._matrix}({self._index+1})'
+        return f'{self._matrix}({self._index+1})'
 
 class ModelMatrix(sym.Matrix,ModelValue):
     """

moyra/symbolic_model.py

@@ -11,6 +11,8 @@
 import pickle
 from .forces import ZeroForce
 from .printing import model as print_model
+from .model_parameters import ModelSymbol, ModelMatrix,ModelMatrixSymbol, ModelValue
+from time import time, ctime
 
 class SymbolicModel:
     """
@@ -283,24 +285,71 @@ def to_matlab_file(self,p,file_dir):
             with open(file_dir+f"{func_name}.m",'w') as file:
                 file.write(self._gen_octave(matrix,p,func_name))
         p.to_matlab_class(file_dir=file_dir)
+
+    def to_matlab_file_linear(self,p,file_dir):
+        mats = self.extract_matrices(p)
+        names = ['A','B','C','D','E']
+        funcs = list(zip([f'get_{i}' for i in names],mats))
+        for func_name,matrix in funcs:
+            with open(file_dir+f"{func_name}.m",'w') as file:
+                file.write(self._gen_octave(matrix,p,func_name))
+        p.to_matlab_class(file_dir=file_dir)
 
     def _gen_octave(self,expr,p,func_name):
         # convert states to non-time dependent variable
         U = sym.Matrix(sym.symbols(f'u_:{p.qs*2}'))
         l = dict(zip(p.x,U))
         expr = me.msubs(expr,l)
 
+        # get parameter replacements
+        param_string = '%% extract required parameters from structure\n\t'
+        matries = []
+        for var in expr.free_symbols:
+            if isinstance(var,ModelValue):
+                if isinstance(var,ModelMatrixSymbol):
+                    if var._matrix not in matries:
+                        param_string += f'{var._matrix} = p.{var._matrix};\n\t'
+                        matries.append(var._matrix)
+                elif isinstance(var,ModelSymbol):
+                    param_string += f'{var.name} = p.{var.name};\n\t'
+                elif isinstance(var,ModelMatrix):
+                    param_string += f'{var._matrix_symbol} = p.{var._matrix_symbol};\n\t'
+
+
+        # split expr into groups
+        replacments, exprs = sym.cse(expr,symbols=(sym.Symbol(f'rep_{i}')for i in range(10000)))
+        if isinstance(expr,tuple):
+            expr = tuple(exprs)
+        elif isinstance(expr,list):
+            expr = exprs
+        else:
+            expr = exprs[0]      
+
+        group_string = '%% create common groups\n\t'
+        for variable, expression in replacments:
+            group_string +=f'{variable} = {sym.printing.octave.octave_code(expression)};\n\t'
+
         # convert to octave string and covert states to vector form
-        out = sym.printing.octave.octave_code(expr)
+        out = '%% create output vector\n\tout = ' + sym.printing.octave.octave_code(expr)
+
+        #convert state vector calls
         my_replace = lambda x: f'U({int(x.group(1))+1})'
         out = re.sub(r"u_(?P<index>\d+)",my_replace,out)
+        group_string = re.sub(r"u_(?P<index>\d+)",my_replace,group_string)
 
         # make the file pretty...
         out = out.replace(',',',...\n\t\t').replace(';',';...\n\t\t')
 
+        file_sig = f'%{func_name.upper()} Auto-generated function from moyra\n\t'
+        file_sig += f'%\n\t'
+        file_sig += f'%\tCreated at : {ctime(time())} \n\t'
+        file_sig += f'%\tCreated with : moyra https://pypi.org/project/moyra/\n\t'
+        file_sig += f'%\n\t'
+
+
         # wrap output in octave function signature
-        signature = f'function out = {func_name}(U,p)'
-        octave_string = signature + '\n\t'+ 'out = ' + out + ';\nend'
+        signature = f'function out = {func_name}(U,p)\n\t'
+        octave_string = signature + file_sig + param_string + group_string + out + ';\nend'
         return octave_string