Merge pull request #299 from ImperialCollegeLondon/dev_hingeaxis

Non-perpendicular hinge axis possible for multibody solver

sharpy/controllers/bladepitchpid.py

@@ -423,7 +423,7 @@ def drive_train_model(self, aero_torque, ini_rot_vel, ini_rot_acc):
 
     def compute_aero_torque(self, beam, struct_tstep):
         # Compute total forces
-        total_forces = np.zeros((6))
+        total_forces = np.zeros(6)
         for ibody in self.settings['blade_num_body']:
             steady, unsteady, grav = struct_tstep.extract_resultants(beam,
                                                       force_type=['steady', 'unsteady', 'grav'],
@@ -436,9 +436,9 @@ def compute_aero_torque(self, beam, struct_tstep):
         hub_node = beam.connectivities[hub_elem, 0]
         hub_pos = struct_tstep.pos[hub_node, :]
 
-        hub_forces = np.zeros((6))
-        hub_forces[0:3] = total_forces[0:3].copy()
-        hub_forces[3:6] = total_forces[3:6] - np.cross(hub_pos, total_forces[0:3])
+        hub_forces = np.zeros(6)
+        hub_forces[:3] = total_forces[:3].copy()
+        hub_forces[3:6] = total_forces[3:6] - np.cross(hub_pos, total_forces[:3])
 
         return hub_forces[5]
 
@@ -450,11 +450,11 @@ def compute_blade_pitch(beam, struct_tstep, tower_ibody=0, blade_ibody=1):
         ielem, inode_in_elem = beam.node_master_elem[tt_node]
         ca0b = algebra.crv2rotation(struct_tstep.psi[ielem, inode_in_elem, :])
         cga0 = algebra.quat2rotation(struct_tstep.mb_quat[tower_ibody, :])
-        zg_tower_top = algebra.multiply_matrices(cga0, ca0b, np.array([0., 0., 1.]))
+        zg_tower_top = cga0 @ ca0b @ np.array([0., 0., 1.])
 
         # blade root
         cga = algebra.quat2rotation(struct_tstep.mb_quat[blade_ibody, :])
-        zg_hub = algebra.multiply_matrices(cga, np.array([0., 0., 1.]))
+        zg_hub = cga @ np.array([0., 0., 1.])
 
         pitch = algebra.angle_between_vectors(zg_tower_top, zg_hub)
         return pitch

sharpy/postproc/aerogridplot.py

@@ -199,7 +199,16 @@ def plot_body(self):
                     counter += 1
                     coords[counter, :] = aero_tstep.zeta[i_surf][:, i_m, i_n]
                     if self.settings['include_rbm']:
-                        coords[counter, :] += struct_tstep.for_pos[0:3]
+                        #TODO: fix for lack of g frame description in nonlineardynamicmultibody.py
+                        if struct_tstep.mb_dict is None:
+                            coords[counter, :] += struct_tstep.for_pos[0:3]
+                        else:
+                            #TODO: uncomment for dynamic trim
+                            # try:
+                            #     cga = algebra.euler2rot([0, self.data.trimmed_values[0], 0])
+                            #     coords[counter, :] += np.dot(cga, struct_tstep.for_pos[0:3])
+                            # except AttributeError:
+                            coords[counter, :] += np.dot(struct_tstep.cga(), struct_tstep.for_pos[0:3])
                     if self.settings['include_forward_motion']:
                         coords[counter, 0] -= self.settings['dt']*self.ts*self.settings['u_inf']
 
@@ -313,7 +322,16 @@ def plot_wake(self):
                     counter += 1
                     coords[counter, :] = self.data.aero.timestep_info[self.ts].zeta_star[i_surf][:, i_m, i_n]
                     if self.settings['include_rbm']:
-                        coords[counter, :] += self.data.structure.timestep_info[self.ts].for_pos[0:3]
+                        #TODO: fix for lack of g frame description in nonlineardynamicmultibody.py
+                        if self.data.structure.timestep_info[self.ts].mb_dict is None:
+                            coords[counter, :] += self.data.structure.timestep_info[self.ts].for_pos[0:3]
+                        else:
+                            #TODO: uncomment for dynamic trim
+                            # try:
+                            #     cga = algebra.euler2rot([0, self.data.trimmed_values[0], 0])
+                            #     coords[counter, :] += np.dot(cga, self.data.structure.timestep_info[self.ts].for_pos[0:3])
+                            # except AttributeError:
+                            coords[counter, :] += np.dot(self.data.structure.timestep_info[self.ts].cga(), self.data.structure.timestep_info[self.ts].for_pos[0:3])                                                        
                     if self.settings['include_forward_motion']:
                         coords[counter, 0] -= self.settings['dt']*self.ts*self.settings['u_inf']
 

sharpy/postproc/beamplot.py

@@ -124,7 +124,6 @@ def write_beam(self, it):
         num_nodes = self.data.structure.num_node
         num_elem = self.data.structure.num_elem
 
-        coords = np.zeros((num_nodes, 3))
         conn = np.zeros((num_elem, 3), dtype=int)
         node_id = np.zeros((num_nodes,), dtype=int)
         elem_id = np.zeros((num_elem,), dtype=int)
@@ -148,6 +147,18 @@ def write_beam(self, it):
         # coordinates of corners
         coords = tstep.glob_pos(include_rbm=self.settings['include_rbm'])
 
+        if tstep.mb_dict is None:
+            pass
+        else:
+            #TODO: fix for lack of g frame description in nonlineardynamicmultibody.py
+            for i_node in range(tstep.num_node):
+                #TODO: uncomment for dynamic trim
+                # try:
+                #     c = algebra.euler2rot([0, self.data.trimmed_values[0], 0])
+                # except AttributeError:
+                c = self.data.structure.timestep_info[0].cga()
+                coords[i_node, :] += np.dot(c, tstep.for_pos[0:3])
+
         # check if I can output gravity forces
         with_gravity = False
         try:
@@ -172,9 +183,6 @@ def write_beam(self, it):
             pass
 
         # count number of arguments
-        postproc_cell_keys = tstep.postproc_cell.keys()
-        postproc_cell_vals = tstep.postproc_cell.values()
-        postproc_cell_scalar = []
         postproc_cell_vector = []
         postproc_cell_6vector = []
         for k, v in tstep.postproc_cell.items():
@@ -189,8 +197,6 @@ def write_beam(self, it):
             else:
                 raise AttributeError('Only scalar and 3-vector types supported in beamplot')
         # count number of arguments
-        postproc_node_keys = tstep.postproc_node.keys()
-        postproc_node_vals = tstep.postproc_node.values()
         postproc_node_scalar = []
         postproc_node_vector = []
         postproc_node_6vector = []

sharpy/postproc/writevariablestime.py

@@ -265,6 +265,28 @@ def write(self, it):
         for ivariable in range(len(self.settings['structure_variables'])):
             if self.settings['structure_variables'][ivariable] == '':
                 continue
+            #TODO: fix for lack of g frame description in nonlineardynamicmultibody.py
+            if tstep.mb_dict is None:
+                var = getattr(tstep, self.settings['structure_variables'][ivariable])
+            else:
+                if self.settings['structure_variables'][ivariable] == 'for_pos':
+                    import sharpy.utils.algebra as ag
+                    #TODO: uncomment for dynamic trim
+                    # try:
+                    #     # import pdb
+                    #     # pdb.set_trace()
+                    #     cga = ag.euler2rot([0, self.data.trimmed_values[0], 0])
+                    #     cag = cga.T
+                    #     tstep.for_pos[0:3] = np.dot(cga,tstep.for_pos[0:3])
+                    #     var = getattr(tstep, self.settings['structure_variables'][ivariable]).copy()
+                    #     tstep.for_pos[0:3] = np.dot(cag,tstep.for_pos[0:3])
+                    # except AttributeError:
+                    t0step = self.data.structure.timestep_info[0]
+                    tstep.for_pos[0:3] = np.dot(t0step.cga(), tstep.for_pos[0:3])
+                    var = getattr(tstep, self.settings['structure_variables'][ivariable]).copy()
+                    tstep.for_pos[0:3] = np.dot(t0step.cag(), tstep.for_pos[0:3])
+                else:
+                    var = getattr(tstep, self.settings['structure_variables'][ivariable])            
             var = getattr(tstep, self.settings['structure_variables'][ivariable])
             num_indices = len(var.shape)
             if num_indices == 1: