# Source code for sharpy.linear.src.lin_utils

```"""
Utilities functions for linear analysis
"""

import numpy as np

# linear uvlm
import sharpy.linear.src.lin_aeroelastic as lin_aeroelastic
import sharpy.linear.src.libss as libss

def comp_tot_force(forces, zeta, zeta_pole=np.zeros((3,))):
""" Compute total force with exact displacements """
Ftot = np.zeros((3,))
Mtot = np.zeros((3,))
for ss in range(len(forces)):
_, Mv, Nv = forces[ss].shape
for mm in range(Mv):
for nn in range(Nv):
arm = zeta[ss][:, mm, nn] - zeta_pole
Mtot += np.cross(arm, forces[ss][:3, mm, nn])
for cc in range(3):
Ftot[cc] += forces[ss][cc, :, :].sum()
return Ftot, Mtot

[docs]class Info():
""" Summarise info about a data point """

def __init__(self, zeta, zeta_dot, u_ext, ftot, mtot, q, qdot,
SSaero=None, SSbeam=None,
Kas=None, Kftot=None, Kmtot=None, Kmtot_disp=None,
self.zeta = zeta
self.zeta_dot = zeta_dot
self.u_ext = u_ext
self.ftot = ftot
self.mtot = mtot
self.q = q
self.qdot = qdot
#
self.SSaero = SSaero
self.SSbeam = SSbeam
self.Kas = Kas
self.Kftot = Kftot
self.Kmtot = Kmtot
self.Kmtot_disp = Kmtot_disp

def solve_linear(Ref, Pert, solve_beam=True):
"""
Given 2 Info() classes associated to a reference linearisation point Ref and a
perturbed state Pert, the method produces in output the prediction at the
Pert state of a linearised model.

The solution is carried on using both the aero and beam input
"""

### define perturbations
dq = Pert.q - Ref.q
dqdot = Pert.qdot - Ref.qdot
dzeta = Pert.zeta - Ref.zeta
dzeta_dot = Pert.zeta_dot - Ref.zeta_dot
du_ext = Pert.u_ext - Ref.u_ext

num_dof_str = len(dq)
dzeta_exp = np.dot(Ref.Kas[:len(dzeta), :num_dof_str], dq)

SSaero = Ref.SSaero
SSbeam = Ref.SSbeam

# zeta in
usta = np.concatenate([dzeta, dzeta_dot, du_ext])

#x_sta = A_steady^-1 dot B u_sta
else:

# y_sta = C x_sta + D u_sta
ysta = np.dot(SSaero.C, xsta) + np.dot(SSaero.D, usta)
ftot_aero = Ref.ftot + np.dot(Ref.Kftot, ysta)
mtot_aero = Ref.mtot + np.dot(Ref.Kmtot, ysta) + np.dot(Ref.Kmtot_disp, dzeta)

#### beam in
if solve_beam:
# warning: we need to add first the contribution due to velocity change!!!
usta_uinf = np.concatenate([0. * dzeta, 0. * dzeta_dot, du_ext])
ysta_uinf = np.dot(SSaero.C, xsta_uinf) + np.dot(SSaero.D, usta_uinf)
usta = np.concatenate([dq, dqdot])
else:
ysta = ysta_uinf + np.dot(SSbeam.C, xsta) + np.dot(SSbeam.D, usta)
ftot_beam = Ref.ftot + np.dot(Ref.Kftot, ysta)
mtot_beam = Ref.mtot + np.dot(Ref.Kmtot, ysta) + np.dot(Ref.Kmtot_disp, dzeta_exp)
else:
ftot_beam, mtot_beam = None, None

return ftot_aero, mtot_aero, ftot_beam, mtot_beam

def extract_from_data(data, assemble=True,
"""
Extract relevant info from data structure. If assemble is True, it will
also generate a linear UVLM and the displacements/velocities gain matrices
"""

### extract aero info - gets info from every panel in every surface
# from an NxM matrix and reshapes it into an K by 1 column vector
tsaero = data.aero.timestep_info[0]
zeta = np.concatenate([tsaero.zeta[ss].reshape(-1, order='C')
for ss in range(tsaero.n_surf)])
zeta_dot = np.concatenate([tsaero.zeta_dot[ss].reshape(-1, order='C')
for ss in range(tsaero.n_surf)])
uext = np.concatenate([tsaero.u_ext[ss].reshape(-1, order='C')
for ss in range(tsaero.n_surf)])
ftot, mtot = comp_tot_force(tsaero.forces, tsaero.zeta, zeta_pole=zeta_pole)

## TEST WHETHER SETTINGS RUN
settings = dict()
settings['LinearUvlm'] = {'dt': 0.1,
'integr_order': 2,
'density': 1.225,
'ScalingDict': {'length': 1.,
'speed': 1.,
'density': 1.}}

### extract structural info
Sol = lin_aeroelastic.LinAeroEla(data, settings)
gebm = Sol.lingebm_str
q = Sol.q
qdot = Sol.dq

### assemble
if assemble is True:
uvlm = Sol.linuvlm
uvlm.assemble_ss()
uvlm.get_total_forces_gain(zeta_pole=zeta_pole)
Sol.get_gebm2uvlm_gains()
Kas = np.block([[Sol.Kdisp, np.zeros((3 * uvlm.Kzeta, gebm.num_dof + 10))],
[Sol.Kvel_disp, Sol.Kvel_vel],
[np.zeros((3 * uvlm.Kzeta, 2 * gebm.num_dof + 20))]])