Repository: ImperialCollegeLondon/sharpy (language: python)

Issue:
Control surface deflection sign convention
Hi all,

This is related to #192 but really is a topic that needs addressing on its own as it is something that we may want to change, although, if we do, will imply backwards changes to all models with symmetric control surfaces.

## TLDR

Control surfaces deflections are not well defined according to a sign convention, which may pose problems in the future for antisymmetric control surfaces (ailerons) and may suppose now a difference in the linear versus nonlinear implementation of the control surfaces. We either change behaviour or continue as is, knowing how it works (a bug becomes a feature if you leave it long enough).

## Background

Control surfaces are implemented in the [aerogrid.py](https://github.com/ImperialCollegeLondon/sharpy/blob/master/sharpy/aero/models/aerogrid.py) file as part of the grid assembly process. 

The UVLM grid is assembled in "strips", where a strip is a chord wise section. These strips are concatenated in a span wise manner resulting in surfaces. 

A strip is generated here: https://github.com/ImperialCollegeLondon/sharpy/blob/41b6cb699fe33167d66cbe0c5ede1bdfaaeb50c4/sharpy/aero/models/aerogrid.py#L482

The whole surface is generated here:
https://github.com/ImperialCollegeLondon/sharpy/blob/41b6cb699fe33167d66cbe0c5ede1bdfaaeb50c4/sharpy/aero/models/aerogrid.py#L205

In the process of generating the strip, control surfaces are done first, in the local B frame by rotating the relevant chord wise section about `xb = [1, 0, 0]`. Then, the section is rotated according to any desired twist or sweep.

## Current process

These strips are concatenated sequentially in increasing node order. However, any rotation on the strip to assemble the surface is performed AFTER the control surface is implemented. Therefore, the control surface is rotated about an axis different from the material frame.

Take, for example the following configuration, where the arrow indicates the local xb, which is aligned with the beam and positive in the increasing node number direction
```
4 -<- 3 -<- 0 ->- 1 ->- 2 
```
The grid strips are generated sequentially for each element. Since the control surfaces are implemented always rotating about the local `[1, 0, 0]` and these are appended, the resulting rotation axis of each CS upon assembly looks like:
```
4 ->- 3 ->- 0 ->- 1 -> 2
```
which is NOT coincident with xb.

## Current behaviour

The current implementation is that, for symmetric control surfaces with a single control surface input (i.e. same `control_surface_id`, control surfaces are deflected in the same direction on both wings, as if they were elevators. Most, if not all, work done in SHARPy has been focused on longitudinal problems so I understand that this is simple as a single input is needed to define an elevator.

### Issues
* The sign convention is not well defined.
* This is not easily implemented in the linear solver, where no strips are generated, and a more clear convention is required.

## Options
1. Leave as is. Backwards compatible and as long as the user is informed on how this works for the nonlinear problem it should be ok. A different convention will be required in the linear solver which may need control surfaces to be defined separately (different `control_surface_id`) as there the convention will be the local `xb`. How would antisymmetric control surfaces (ailerons) be implemented with the current implementation and a single input is yet to be determined.
2. Change the behaviour, this can be done by moving where the control surface is implemented in https://github.com/ImperialCollegeLondon/sharpy/blob/41b6cb699fe33167d66cbe0c5ede1bdfaaeb50c4/sharpy/aero/models/aerogrid.py#L482
moving the `#control surface implementation` after `#transformation to beam and beam prime`.  The result is 
![image](https://user-images.githubusercontent.com/35684407/165939685-b98e229e-39b7-4429-a72a-0f71e435539b.png)
If this is done, this will affect all previous models as effectively all previously implemented control surfaces will deflect anti-symmetrically. The current behaviour (i.e. elevator-like) could be retained by modifying the models `xb` as 
```
3 ->- 4 ->- 0 ->- 1 -> 2
```
The advantage of this being that the sign convention is clearly defined, and that the same approach is retained in the linear solver.
3. A final option would be to leave as is, introducing some kind of "linkage" between control surfaces. This could mean a -1 gain from one to another to simulate ailerons in the current nonlinear implementation.

## Please comment :)
A bit of a longwinded explanation but I think that should be tackled as the code goes forward. Probably the easiest short term option is 1) (and I'll see what to do for the linear solver in relation to #192) but we may want to look towards option 2) for future releases.


Discussion / hints:
Brilliant explanation. A few thoughts after quick look to the code:

1. The part of the code that makes the trick for the symmetric behaviour of the control surface on right and left wings is the definition of this rot_angle https://github.com/ImperialCollegeLondon/sharpy/blob/41b6cb699fe33167d66cbe0c5ede1bdfaaeb50c4/sharpy/aero/models/aerogrid.py#L557
this is also what makes definition of twist consistent. 
I would add another index to the variable `control_surface_type`, -1 for example, for antisymmetric behaviour of ailerons, and depending on that apply the Crot before or after the control surface rotation. Thus maintaining backwards compatibility but having this added feature in the future (a table in the docs with the options would suffice for establishing a convention for users). I would not change xb vectors which are well defined towards increasing node-index and are sort of lower level. For the linear, we will need to have a chat. 

2. We would need to check whether this Crot applies to complex configurations with structural sweep, dihedral, etc.  It is part of a bigger discussion regarding the symmetric modelling of the aircraft since some issues are appearing when inputting fully populated stiffness and mass matrices on these models. It might be worth opening a github project for this, I will arrange something when freer. 

3. Before doing this we need to merge Arturo's last commit where he put the option to apply sweep before twist...

4. @kccwing, this is a very relevant discussion for your project regarding the semi-aeroelastic hinge modelling. When you sort out the multibody part, the next challenge is to have an aero surface that is continuous for the different directions of the hinge... it could entail some additions to this aerogrid.py

Yep, the thing is that control surfaces are applied currently before the definition of `rot_angle`, thus they are applied using `rotation_3dx` in which the `x` axis used for the rotation is not necessarily `xb`. I wouldn't go into changing `xb`, but only moving where the CS is applied.

Practically, most applications benefit from the symmetric CS behaviour as is, therefore maybe the best option is to add the extra input to tell the code whether the CS deflects anti-symmetrically.

Agree on 2). I had experience with symmetric modelling with non-diagonal stiffness matrices and it was tricky.

Also agree on 3). Although the control surface behaviour is not affected by it

Fix patch:
diff --git a/sharpy/aero/models/aerogrid.py b/sharpy/aero/models/aerogrid.py
index 0030b56b9..399ce6c0c 100644
--- a/sharpy/aero/models/aerogrid.py
+++ b/sharpy/aero/models/aerogrid.py
@@ -188,13 +188,11 @@ def calculate_dimensions(self):
                     nodes_in_surface[i_surf].append(i_global_node)
                 if self.aero_dict['aero_node'][i_global_node]:
                     self.aero_dimensions[i_surf, 1] += 1
-
         # accounting for N+1 nodes -> N panels
-        self.aero_dimensions[:, 1] -= 1
-
+        self.aero_dimensions[:,1] -= 1
+        
         self.aero_dimensions_star = self.aero_dimensions.copy()
-        for i_surf in range(self.n_surf):
-            self.aero_dimensions_star[i_surf, 0] = self.aero_settings['mstar']
+        self.aero_dimensions_star[:, 0] = self.aero_settings['mstar']
 
     def add_timestep(self):
         try:
diff --git a/sharpy/aero/utils/airfoilpolars.py b/sharpy/aero/utils/airfoilpolars.py
index 40d55f707..514faadd1 100644
--- a/sharpy/aero/utils/airfoilpolars.py
+++ b/sharpy/aero/utils/airfoilpolars.py
@@ -1,6 +1,7 @@
 import numpy as np
 # import sharpy.utils.algebra as algebra
 from sharpy.utils.constants import deg2rad
+from scipy.interpolate import interp1d
 
 
 class Polar:
@@ -23,7 +24,7 @@ def initialise(self, table):
         """
         # Store the table
         if (np.diff(table[:, 0]) > 0.).all():
-            self.table = table
+            self.table = table[~np.isnan(table).any(axis=1), :]
         else:
             raise RuntimeError("ERROR: angles of attack not ordered")
 
@@ -51,11 +52,15 @@ def initialise(self, table):
                 iaoacl0 = imin
         self.aoa_cl0_deg = matches[iaoacl0]
 
+        self.cl_interp = interp1d(self.table[:, 0], self.table[:, 1])
+        self.cd_interp = interp1d(self.table[:, 0], self.table[:, 2])
+        self.cm_interp = interp1d(self.table[:, 0], self.table[:, 3])
+
     def get_coefs(self, aoa_deg):
 
-        cl = np.interp(aoa_deg, self.table[:, 0], self.table[:, 1])
-        cd = np.interp(aoa_deg, self.table[:, 0], self.table[:, 2])
-        cm = np.interp(aoa_deg, self.table[:, 0], self.table[:, 3])
+        cl = self.cl_interp(aoa_deg)
+        cd = self.cd_interp(aoa_deg)
+        cm = self.cm_interp(aoa_deg)
 
         return cl, cd, cm
 
@@ -63,6 +68,7 @@ def get_aoa_deg_from_cl_2pi(self, cl):
 
         return cl/2/np.pi/deg2rad + self.aoa_cl0_deg
 
+
     def redefine_aoa(self, new_aoa):
 
         naoa = len(new_aoa)
diff --git a/sharpy/aero/utils/utils.py b/sharpy/aero/utils/utils.py
index 7a2547495..93757372f 100644
--- a/sharpy/aero/utils/utils.py
+++ b/sharpy/aero/utils/utils.py
@@ -130,38 +130,27 @@ def span_chord(i_node_surf, zeta):
     return dir_span, span, dir_chord, chord
 
 
-def find_aerodynamic_solver(settings):
+def find_aerodynamic_solver_settings(settings):
     """
-    Retrieves the name and settings of the first aerodynamic solver used in the solution ``flow``.
+    Retrieves the settings of the first aerodynamic solver used in the solution ``flow``.
+
 
     Args:
         settings (dict): SHARPy settings (usually found in ``data.settings`` )
 
     Returns:
-        tuple: Aerodynamic solver name and solver settings
+        tuple: Aerodynamic solver settings
     """
     flow = settings['SHARPy']['flow']
-    # Look for the aerodynamic solver
-    if 'StaticUvlm' in flow:
-        aero_solver_name = 'StaticUvlm'
-        aero_solver_settings = settings['StaticUvlm']
-    elif 'StaticCoupled' in flow:
-        aero_solver_name = settings['StaticCoupled']['aero_solver']
-        aero_solver_settings = settings['StaticCoupled']['aero_solver_settings']
-    elif 'StaticCoupledRBM' in flow:
-        aero_solver_name = settings['StaticCoupledRBM']['aero_solver']
-        aero_solver_settings = settings['StaticCoupledRBM']['aero_solver_settings']
-    elif 'DynamicCoupled' in flow:
-        aero_solver_name = settings['DynamicCoupled']['aero_solver']
-        aero_solver_settings = settings['DynamicCoupled']['aero_solver_settings']
-    elif 'StepUvlm' in flow:
-        aero_solver_name = 'StepUvlm'
-        aero_solver_settings = settings['StepUvlm']
-    else:
-        raise KeyError("ERROR: aerodynamic solver not found")
-
-    return aero_solver_name, aero_solver_settings
+    for solver_name in ['StaticUvlm', 'StaticCoupled', 'DynamicCoupled', 'StepUvlm']:
+        if solver_name in flow:
+            aero_solver_settings = settings[solver_name]
+            if 'aero_solver' in settings[solver_name].keys():
+                aero_solver_settings = aero_solver_settings['aero_solver_settings']
+                
+            return aero_solver_settings
 
+    raise KeyError("ERROR: Aerodynamic solver not found.")
 
 def find_velocity_generator(settings):
     """
@@ -174,7 +163,7 @@ def find_velocity_generator(settings):
     Returns:
         tuple: velocity generator name and velocity generator settings
     """
-    aero_solver_name, aero_solver_settings = find_aerodynamic_solver(settings)
+    aero_solver_settings = find_aerodynamic_solver_settings(settings)
 
     vel_gen_name = aero_solver_settings['velocity_field_generator']
     vel_gen_settings = aero_solver_settings['velocity_field_input']
diff --git a/sharpy/generators/gustvelocityfield.py b/sharpy/generators/gustvelocityfield.py
index 2e8b579e7..e13157a53 100644
--- a/sharpy/generators/gustvelocityfield.py
+++ b/sharpy/generators/gustvelocityfield.py
@@ -19,6 +19,7 @@
 from abc import ABCMeta
 import sharpy.utils.generator_interface as generator_interface
 import sharpy.utils.settings as settings
+from scipy.interpolate import interp1d
 
 dict_of_gusts = {}
 
@@ -88,6 +89,10 @@ class one_minus_cos(BaseGust):
     settings_default['gust_intensity'] = 0.0
     settings_description['gust_intensity'] = 'Intensity of the gust :math:`u_{de}`.'
 
+    settings_types['gust_component'] = 'int'
+    settings_default['gust_component'] = 2
+    settings_description['gust_component'] = 'Component of the gust velocity in the G-frame (x,y,z)->(0,1,2).'
+
     setting_table = settings.SettingsTable()
     __doc__ += setting_table.generate(settings_types, settings_default, settings_description,
                                       header_line=doc_settings_description)
@@ -104,7 +109,7 @@ def gust_shape(self, x, y, z, time=0):
         if x > 0.0 or x < -gust_length:
             return vel
 
-        vel[2] = (1.0 - np.cos(2.0 * np.pi * x / gust_length)) * gust_intensity * 0.5
+        vel[self.settings['gust_component']] = (1.0 - np.cos(2.0 * np.pi * x / gust_length)) * gust_intensity * 0.5
         return vel
 
 
@@ -135,6 +140,10 @@ class DARPA(BaseGust):
     settings_default['span'] = 0.
     settings_description['span'] = 'Wing span'
 
+    settings_types['gust_component'] = 'int'
+    settings_default['gust_component'] = 2
+    settings_description['gust_component'] = 'Component of the gust velocity in the G-frame (x,y,z)->(0,1,2).'
+
     setting_table = settings.SettingsTable()
     __doc__ += setting_table.generate(settings_types, settings_default, settings_description,
                                       header_line=doc_settings_description)
@@ -152,8 +161,8 @@ def gust_shape(self, x, y, z, time=0):
         if x > 0.0 or x < -gust_length:
             return vel
 
-        vel[2] = (1.0 - np.cos(2.0 * np.pi * x / gust_length)) * gust_intensity * 0.5
-        vel[2] *= -np.cos(y / span * np.pi)
+        vel[self.settings['gust_component']] = (1.0 - np.cos(2.0 * np.pi * x / gust_length)) * gust_intensity * 0.5
+        vel[self.settings['gust_component']] *= -np.cos(y / span * np.pi)
         return vel
 
 
@@ -180,44 +189,9 @@ class continuous_sin(BaseGust):
     settings_default['gust_intensity'] = 0.0
     settings_description['gust_intensity'] = 'Intensity of the gust'
 
-    setting_table = settings.SettingsTable()
-    __doc__ += setting_table.generate(settings_types, settings_default, settings_description,
-                                      header_line=doc_settings_description)
-
-    def initialise(self, in_dict, restart=False):
-        self.settings = in_dict
-        settings.to_custom_types(self.settings, self.settings_types, self.settings_default)
-
-    def gust_shape(self, x, y, z, time=0):
-        gust_length = self.settings['gust_length']
-        gust_intensity = self.settings['gust_intensity']
-
-        vel = np.zeros((3,))
-        if x > 0.0:
-            return vel
-
-        vel[2] = 0.5 * gust_intensity * np.sin(2 * np.pi * x / gust_length)
-        return vel
-
-
-@gust
-class lateral_one_minus_cos(BaseGust):
-    r"""
-    This gust can be used by using the setting ``gust_shape = 'lateral 1-cos'`` in :class:`GustVelocityField`.
-    """
-    gust_id = 'lateral 1-cos'
-
-    settings_types = dict()
-    settings_default = dict()
-    settings_description = dict()
-
-    settings_types['gust_length'] = 'float'
-    settings_default['gust_length'] = 0.0
-    settings_description['gust_length'] = 'Length of gust'
-
-    settings_types['gust_intensity'] = 'float'
-    settings_default['gust_intensity'] = 0.0
-    settings_description['gust_intensity'] = 'Intensity of the gust'
+    settings_types['gust_component'] = 'int'
+    settings_default['gust_component'] = 2
+    settings_description['gust_component'] = 'Component of the gust velocity in the G-frame (x,y,z)->(0,1,2).'
 
     setting_table = settings.SettingsTable()
     __doc__ += setting_table.generate(settings_types, settings_default, settings_description,
@@ -232,23 +206,22 @@ def gust_shape(self, x, y, z, time=0):
         gust_intensity = self.settings['gust_intensity']
 
         vel = np.zeros((3,))
-        if x > 0.0 or x < -gust_length:
+        if x > 0.0:
             return vel
 
-        vel[1] = (1.0 - np.cos(2.0 * np.pi * x / gust_length)) * gust_intensity * 0.5
+        vel[self.settings['gust_component']] = 0.5 * gust_intensity * np.sin(2 * np.pi * x / gust_length)
         return vel
 
 
 @gust
-class time_varying(BaseGust):
+class time_varying_global(BaseGust):
     r"""
-    The inflow velocity changes with time but it is uniform in space. It is read from a 4 column file:
-
-    .. math:: time[s] \Delta U_x \Delta U_y \Delta U_z
+    Similar to the previous one but the velocity changes instanteneously in the whole flow field. It is not fed
+    into the solid.
 
-    This gust can be used by using the setting ``gust_shape = 'time varying'`` in :class:.`GustVelocityField`.
+    This gust can be used by using the setting ``gust_shape = 'time varying global'`` in :class:`GustVelocityField`.
     """
-    gust_id = 'time varying'
+    gust_id = 'time varying global'
 
     settings_types = dict()
     settings_default = dict()
@@ -256,7 +229,11 @@ class time_varying(BaseGust):
 
     settings_types['file'] = 'str'
     settings_default['file'] = ''
-    settings_description['file'] = 'File with the information'
+    settings_description['file'] = 'File with the information (only for time varying)'
+
+    settings_types['gust_component'] = ['list(int)', 'int']
+    settings_default['gust_component'] = [0, 1, 2]
+    settings_description['gust_component'] = 'Component of the gust velocity in the G-frame to be considered (x,y,z)->(0,1,2).'
 
     setting_table = settings.SettingsTable()
     __doc__ += setting_table.generate(settings_types, settings_default, settings_description,
@@ -264,69 +241,52 @@ class time_varying(BaseGust):
 
     def __init__(self):
         super().__init__()
-
         self.file_info = None
+        self.list_interpolated_velocity_field_functions = []
 
     def initialise(self, in_dict, restart=False):
         self.settings = in_dict
         settings.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         self.file_info = np.loadtxt(self.settings['file'])
+        self.initialise_interpolation_functions()
+    
+    def initialise_interpolation_functions(self):
+        for idim in self.settings['gust_component']:
+            self.list_interpolated_velocity_field_functions.append(interp1d(self.file_info[:, 0], self.file_info[:, idim+1], 
+                                                                            bounds_error=False,fill_value="extrapolate"))
 
     def gust_shape(self, x, y, z, time=0):
         vel = np.zeros((3,))
-        d = np.dot(np.array([x, y, z]), self.u_inf_direction)
-        if d > 0.0:
-            return vel
-
-        vel[0] = np.interp(d, -self.file_info[::-1, 0] * self.u_inf, self.file_info[::-1, 1])
-        vel[1] = np.interp(d, -self.file_info[::-1, 0] * self.u_inf, self.file_info[::-1, 2])
-        vel[2] = np.interp(d, -self.file_info[::-1, 0] * self.u_inf, self.file_info[::-1, 3])
+        for counter, idim in enumerate(self.settings['gust_component']):
+            vel[idim] = self.list_interpolated_velocity_field_functions[counter](time)
         return vel
 
 
 @gust
-class time_varying_global(BaseGust):
+class time_varying(time_varying_global):
     r"""
-    Similar to the previous one but the velocity changes instanteneously in the whole flow field. It is not fed
-    into the solid.
-
-    This gust can be used by using the setting ``gust_shape = 'time varying global'`` in :class:`GustVelocityField`.
-    """
-    gust_id = 'time varying global'
-
-    settings_types = dict()
-    settings_default = dict()
-    settings_description = dict()
-
-    settings_types['file'] = 'str'
-    settings_default['file'] = ''
-    settings_description['file'] = 'File with the information (only for time varying)'
-
-    setting_table = settings.SettingsTable()
-    __doc__ += setting_table.generate(settings_types, settings_default, settings_description,
-                                      header_line=doc_settings_description)
-
-    def __init__(self):
-        super().__init__()
-
-        self.file_info = None
+    The inflow velocity changes with time but it is uniform in space. It is read from a 4 column file:
 
-    def initialise(self, in_dict, restart=False):
-        self.settings = in_dict
-        settings.to_custom_types(self.settings, self.settings_types, self.settings_default)
+    .. math:: time[s] \Delta U_x \Delta U_y \Delta U_z
 
-        self.file_info = np.loadtxt(self.settings['file'])
+    This gust can be used by using the setting ``gust_shape = 'time varying'`` in :class:.`GustVelocityField`.
+    """
+    gust_id = 'time varying'
+    
+    def initialise_interpolation_functions(self):
+        for idim in self.settings['gust_component']:
+            self.list_interpolated_velocity_field_functions.append(interp1d(-self.file_info[::-1, 0] * self.u_inf, self.file_info[::-1, idim+1], 
+                                                                            bounds_error=False,fill_value="extrapolate"))
 
     def gust_shape(self, x, y, z, time=0):
         vel = np.zeros((3,))
-
-        vel[0] = np.interp(time, self.file_info[:, 0], self.file_info[:, 1])
-        vel[1] = np.interp(time, self.file_info[:, 0], self.file_info[:, 2])
-        vel[2] = np.interp(time, self.file_info[:, 0], self.file_info[:, 3])
+        d = np.dot(np.array([x, y, z]), self.u_inf_direction)
+        if d <= 0.0:       
+            for counter, idim in enumerate(self.settings['gust_component']):
+                vel[idim] = self.list_interpolated_velocity_field_functions[counter](d)
         return vel
-
-
+       
 @gust
 class span_sine(BaseGust):
     r"""
diff --git a/sharpy/generators/polaraeroforces.py b/sharpy/generators/polaraeroforces.py
index c68d78880..b17e59ff0 100644
--- a/sharpy/generators/polaraeroforces.py
+++ b/sharpy/generators/polaraeroforces.py
@@ -1,4 +1,5 @@
 import numpy as np
+import os
 import sharpy.utils.generator_interface as generator_interface
 import sharpy.utils.settings as settings
 import sharpy.utils.algebra as algebra
@@ -69,7 +70,7 @@ class PolarCorrection(generator_interface.BaseGenerator):
                                                 'simply from polar derived data, i.e. the polars Cm and the moments' \
                                                 'arising from the lift and drag (derived from the polar) contribution. ' \
                                                 'Else, it will add the polar Cm to the moment already computed by ' \
-                                                'SHARPy.'
+                                                'SHARPy.'    
 
     settings_types['add_rotation'] = 'bool'
     settings_default['add_rotation'] = False
@@ -83,6 +84,19 @@ class PolarCorrection(generator_interface.BaseGenerator):
     settings_default['centre_rot_g'] = [0., 0., 0.] 
     settings_description['centre_rot_g'] = 'Centre of rotation in G FoR. Only used if add_rotation = True'
 
+    settings_types['skip_surfaces'] = 'list(int)'
+    settings_default['skip_surfaces'] = []
+    settings_description['skip_surfaces'] = 'Surfaces on which force correction is skipped.'
+
+
+    settings_types['aoa_cl0'] = 'list(float)'
+    settings_default['aoa_cl0'] = []
+    settings_description['aoa_cl0'] = 'Angle of attack for which zero lift is achieved specified in deg for each airfoil.'
+    
+    settings_types['write_induced_aoa'] = 'bool'
+    settings_default['write_induced_aoa'] = False
+    settings_description['write_induced_aoa'] = 'Write induced aoa of each node to txt file.'
+
     settings_table = settings.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options,
                                        header_line='This generator takes in the following settings.')
@@ -94,6 +108,8 @@ def __init__(self):
         self.structure = None
         self.rho = None
         self.vortex_radius = None
+        self.n_node = None
+        self.flag_node_shared_by_multiple_surfaces = None
 
     def initialise(self, in_dict, **kwargs):
         self.settings = in_dict
@@ -101,8 +117,19 @@ def initialise(self, in_dict, **kwargs):
 
         self.aero = kwargs.get('aero')
         self.structure = kwargs.get('structure')
+        self.n_node = self.structure.num_node
         self.rho = kwargs.get('rho')
         self.vortex_radius = kwargs.get('vortex_radius', 1e-6)
+        self.list_aoa_cl0 = self.settings['aoa_cl0']
+        self.cd_from_cl = self.settings['cd_from_cl']
+        self.folder = kwargs.get('output_folder') + '/aoa_induced/'
+
+        if not self.cd_from_cl and len(self.list_aoa_cl0) == 0:
+            # compute aoa for cl0 if not specified in settings
+            self.compute_aoa_cl0_from_airfoil_data(self.aero)
+
+        self.check_for_special_cases(self.aero)
+
 
     def generate(self, **params):
         """
@@ -118,13 +145,15 @@ def generate(self, **params):
         aero_kstep = params['aero_kstep']
         structural_kstep = params['structural_kstep']
         struct_forces = params['struct_forces']
+        ts = params['ts']
 
         aerogrid = self.aero
         structure = self.structure
         rho = self.rho
         correct_lift = self.settings['correct_lift']
-        cd_from_cl = self.settings['cd_from_cl']
         moment_from_polar = self.settings['moment_from_polar']
+        
+        list_aoa_induced = []
 
         aero_dict = aerogrid.aero_dict
         if aerogrid.polars is None:
@@ -143,67 +172,70 @@ def generate(self, **params):
                 ielem, inode_in_elem = structure.node_master_elem[inode]
                 iairfoil = aero_dict['airfoil_distribution'][ielem, inode_in_elem]
                 isurf = aerogrid.struct2aero_mapping[inode][0]['i_surf']
-                i_n = aerogrid.struct2aero_mapping[inode][0]['i_n']
-                N = aerogrid.aero_dimensions[isurf, 1]
-                polar = aerogrid.polars[iairfoil]
-                cab = algebra.crv2rotation(structural_kstep.psi[ielem, inode_in_elem, :])
-                cgb = np.dot(cga, cab)
-
-                if not cd_from_cl:
-                    airfoil_coords = aerogrid.aero_dict['airfoils'][str(aero_dict['airfoil_distribution'][ielem, inode_in_elem])]
-
-                dir_span, span, dir_chord, chord = span_chord(i_n, aero_kstep.zeta[isurf])
-
-                # Define the relative velocity and its direction
-                urel, dir_urel = magnitude_and_direction_of_relative_velocity(structural_kstep.pos[inode, :],
-                                                                              structural_kstep.pos_dot[inode, :],
-                                                                              structural_kstep.for_vel[:],
-                                                                              cga,
-                                                                              aero_kstep.u_ext[isurf][:, :, i_n],
-                                                                              self.settings['add_rotation'],
-                                                                              self.settings['rot_vel_g'],
-                                                                              self.settings['centre_rot_g'],)
-
-                # Coefficient to change from aerodynamic coefficients to forces (and viceversa)
-                coef = 0.5 * rho * np.linalg.norm(urel) ** 2 * chord * span
-
-                # Stability axes - projects forces in B onto S
-                c_bs = local_stability_axes(cgb.T.dot(dir_urel), cgb.T.dot(dir_chord))
-                forces_s = c_bs.T.dot(struct_forces[inode, :3])
-                moment_s = c_bs.T.dot(struct_forces[inode, 3:])
-                drag_force = forces_s[0]
-                lift_force = forces_s[2]
-
-                # Compute the associated lift
-                cl = np.sign(lift_force) * np.linalg.norm(lift_force) / coef
-                cd_sharpy = np.linalg.norm(drag_force) / coef
-
-                if cd_from_cl:
-                    # Compute the drag from the UVLM computed lift
-                    cd, cm = polar.get_cdcm_from_cl(cl)
-
-                else:
-                    # Compute L, D, M from polar depending on:
-                    # ii) Compute the effective angle of attack from potential flow theory. The local lift curve
-                    # slope is 2pi and the zero-lift angle of attack is given by thin airfoil theory. From this,
-                    # the effective angle of attack is computed for the section and includes 3D effects.
-                    aoa_0cl = get_aoacl0_from_camber(airfoil_coords[:, 0], airfoil_coords[:, 1])
-                    aoa = cl / 2 / np.pi + aoa_0cl
-                    # Compute the coefficients associated to that angle of attack
-                    cl_polar, cd, cm = polar.get_coefs(aoa)
-
-                    if correct_lift:
-                        # Use polar generated CL rather than UVLM computed CL
-                        cl = cl_polar
-
-                # Recompute the forces based on the coefficients (side force is uncorrected)
-                forces_s[0] += cd * coef  # add viscous drag to induced drag from UVLM
-                forces_s[2] = cl * coef
-
-                new_struct_forces[inode, 0:3] = c_bs.dot(forces_s)
-
-                # Pitching moment
-                if moment_from_polar:
+                if isurf not in self.settings['skip_surfaces']:
+                    i_n = aerogrid.struct2aero_mapping[inode][0]['i_n']
+                    polar = aerogrid.polars[iairfoil]
+                    cab = algebra.crv2rotation(structural_kstep.psi[ielem, inode_in_elem, :])
+                    cgb = np.dot(cga, cab)
+
+                    if not self.cd_from_cl:
+                        airfoil = str(aero_dict['airfoil_distribution'][ielem, inode_in_elem])
+                        aoa_0cl = self.list_aoa_cl0[int(airfoil)]
+
+                    # computing surface area of panels contributing to force
+                    dir_span, span, dir_chord, chord = span_chord(i_n, aero_kstep.zeta[isurf])
+                    area = span * chord
+                    area = self.correct_surface_area(inode, aerogrid.struct2aero_mapping, aero_kstep.zeta, area)
+               
+                    # Define the relative velocity and its direction
+                    urel, dir_urel = magnitude_and_direction_of_relative_velocity(structural_kstep.pos[inode, :],
+                                                                                structural_kstep.pos_dot[inode, :],
+                                                                                structural_kstep.for_vel[:],
+                                                                                cga,
+                                                                                aero_kstep.u_ext[isurf][:, :, i_n],
+                                                                                self.settings['add_rotation'],
+                                                                                self.settings['rot_vel_g'],
+                                                                                self.settings['centre_rot_g'],)
+                    # Coefficient to change from aerodynamic coefficients to forces (and viceversa)
+                    coef = 0.5 * rho * np.linalg.norm(urel) ** 2 * area
+                    # Stability axes - projects forces in B onto S
+                    c_bs = local_stability_axes(cgb.T.dot(dir_urel), cgb.T.dot(dir_chord))
+                    forces_s = c_bs.T.dot(struct_forces[inode, :3])
+                    moment_s = c_bs.T.dot(struct_forces[inode, 3:])
+                    drag_force = forces_s[0]
+                    lift_force = forces_s[2]
+                    # Compute the associated lift
+                    cl = np.sign(lift_force) * np.linalg.norm(lift_force) / coef
+                    cd_sharpy = np.linalg.norm(drag_force) / coef
+
+                    if self.cd_from_cl:
+                        # Compute the drag from the UVLM computed lift
+                        cd, cm = polar.get_cdcm_from_cl(cl)
+
+                    else:
+                        """
+                        Compute L, D, M from polar depending on:
+                        ii) Compute the effective angle of attack from potential flow theory or specified it as setting
+                        input. The local lift curve slope is 2pi and the zero-lift angle of attack is given by thin 
+                        airfoil theory or specified it as setting input. From this, the effective angle of attack is 
+                        computed for the section and includes 3D effects.
+                        """
+                        aoa = cl / 2 / np.pi + aoa_0cl 
+                        list_aoa_induced.append(aoa)
+                        # Compute the coefficients associated to that angle of attack
+                        cl_polar, cd, cm = polar.get_coefs(aoa)
+                        
+                        if correct_lift:
+                            # Use polar generated CL rather than UVLM computed CL
+                            cl = cl_polar
+
+                    # Recompute the forces based on the coefficients (side force is uncorrected)
+                    forces_s[0] += cd * coef  # add viscous drag to induced drag from UVLM
+                    forces_s[2] = cl * coef
+
+                    new_struct_forces[inode, 0:3] = c_bs.dot(forces_s)
+
+                    # Pitching moment
                     # The panels are shifted by 0.25 of a panel aft from the leading edge
                     panel_shift = 0.25 * (aero_kstep.zeta[isurf][:, 1, i_n] - aero_kstep.zeta[isurf][:, 0, i_n])
                     ref_point = aero_kstep.zeta[isurf][:, 0, i_n] + 0.25 * chord * dir_chord - panel_shift
@@ -216,6 +248,22 @@ def generate(self, **params):
                     moment_polar_drag = algebra.cross3(c_bs.T.dot(arm), cd * dir_urel * coef)  # in S frame
                     moment_s += moment_polar_drag
 
+                    # Pitching moment
+                    if moment_from_polar:
+                        # The panels are shifted by 0.25 of a panel aft from the leading edge
+                        panel_shift = 0.25 * (aero_kstep.zeta[isurf][:, 1, i_n] - aero_kstep.zeta[isurf][:, 0, i_n])
+                        ref_point = aero_kstep.zeta[isurf][:, 0, i_n] + 0.25 * chord * dir_chord - panel_shift
+                        new_struct_forces[inode, 3:6] = c_bs.dot(moment_s)
+
+                        # viscous contribution (pure moment)
+                        moment_s[1] += cm * coef * chord
+
+
+                        # moment due to drag
+                        arm = cgb.T.dot(ref_point - pos_g[inode])  # in B frame
+                        moment_polar_drag = algebra.cross3(c_bs.T.dot(arm), cd * dir_urel * coef)  # in S frame
+                        moment_s += moment_polar_drag
+
                     # moment due to lift (if corrected)
                     if correct_lift and moment_from_polar:
                         # add moment from scratch: cm_polar + cm_drag_polar + cl_lift_polar
@@ -227,8 +275,85 @@ def generate(self, **params):
 
                     new_struct_forces[inode, 3:6] = c_bs.dot(moment_s)
 
+        if self.settings['write_induced_aoa']:
+            self.write_induced_aoa_of_each_node(ts, list_aoa_induced)
+
         return new_struct_forces
+    
+
+    def correct_surface_area(self, inode, struct2aero_mapping, zeta_ts, area):
+        '''
+        Corrects the surface area if the structural node is shared  by multiple surfaces. 
+
+        For example, when the wing is split into right and left wing both surfaces share the center node.
+        Necessary for cl calculation as the force on the node is already the sum of the forces generated 
+        at the adjacent panels of each surface.
+        
+        Args:
+            inode (int): global node id
+            struct2aero_mapping (list of dicts): maps the structural (global) nodes to aero surfaces and nodes
+            zeta_ts (array): zeta of current aero timestep
 
+        Returns:
+            float: corrected surface area of other surfaces
+        '''
+        if self.flag_shared_node_by_surfaces[inode]:
+            n_surfaces_shared_by_node = len(struct2aero_mapping[inode])
+            # add area for all other surfaces connected to this node
+            for isurf in range(1,n_surfaces_shared_by_node):
+                shared_surf = struct2aero_mapping[inode][isurf]['i_surf']
+                i_n_shared_surf = struct2aero_mapping[inode][isurf]['i_n']
+                _, span_shared_surf, _, chord_shared_surf = span_chord(i_n_shared_surf, zeta_ts[shared_surf])
+                area += span_shared_surf * chord_shared_surf
+        return area
+        
+    def check_for_special_cases(self, aerogrid):
+        '''
+        Checks if the outboard node is shared by multiple surfaces. 
+
+        Args:
+            aerogrid :class:`~sharpy.aero.models.AerogridLoader
+        '''
+        # check if outboard node of aerosurface
+        self.flag_shared_node_by_surfaces = np.zeros((self.n_node,1))
+        for inode in range(self.n_node):
+            if aerogrid.aero_dict['aero_node'][inode]:
+                i_n = aerogrid.struct2aero_mapping[inode][0]['i_n']                
+                isurf = aerogrid.struct2aero_mapping[inode][0]['i_surf']
+                N = aerogrid.aero_dimensions[isurf, 1]
+                if i_n in [0, N]:
+                    if len(aerogrid.struct2aero_mapping[inode]) > 1:
+                        self.flag_shared_node_by_surfaces[inode] = 1
+      
+
+
+    def write_induced_aoa_of_each_node(self,ts, list_aoa_induced):
+        '''
+        Writes induced aoa of each node to txt file for each timestep. 
+
+        Args:
+            ts (int): simulation timestep 
+            list_aoa_induced (list(float)): list with induced aoa of each node
+        '''
+        if ts == 0 and not os.path.exists(self.folder):
+            os.makedirs(self.folder)
+            
+        np.savetxt(self.folder + '/aoa_induced_ts_{}.txt'.format(ts),
+                   np.transpose(np.transpose(np.array(list_aoa_induced))),
+                   fmt='%10e',
+                   delimiter=',',
+                   header='aoa_induced',
+                   comments='#')
+                   
+
+    def compute_aoa_cl0_from_airfoil_data(self, aerogrid):
+        """
+         Computes the angle of attack for which zero lift is achieved for every airfoil
+        """
+        self.list_aoa_cl0 = np.zeros((len(aerogrid.aero_dict['airfoils']),1))
+        for i, airfoil in enumerate(aerogrid.aero_dict['airfoils']):
+            airfoil_coords = aerogrid.aero_dict['airfoils'][airfoil]
+            self.list_aoa_cl0[i] = get_aoacl0_from_camber(airfoil_coords[:, 0], airfoil_coords[:, 1])
 
 @generator_interface.generator
 class EfficiencyCorrection(generator_interface.BaseGenerator):
diff --git a/sharpy/generators/straightwake.py b/sharpy/generators/straightwake.py
index 095f76a30..2bd31af85 100644
--- a/sharpy/generators/straightwake.py
+++ b/sharpy/generators/straightwake.py
@@ -41,19 +41,19 @@ class StraightWake(generator_interface.BaseGenerator):
     settings_default['dt'] = 0.1 # None
     settings_description['dt'] = 'Time step'
 
-    settings_types['dx1'] = 'float'
+    settings_types['dx1'] =  ['float', 'list(float)']
     settings_default['dx1'] = -1.0
     settings_description['dx1'] = 'Size of the first wake panel'
 
-    settings_types['ndx1'] = 'int'
+    settings_types['ndx1'] =  ['int', 'list(int)']
     settings_default['ndx1'] = 1
     settings_description['ndx1'] = 'Number of panels with size ``dx1``'
 
-    settings_types['r'] = 'float'
+    settings_types['r'] = ['float', 'list(float)']
     settings_default['r'] = 1.
     settings_description['r'] = 'Growth rate after ``ndx1`` panels'
 
-    settings_types['dxmax'] = 'float'
+    settings_types['dxmax'] =  ['float', 'list(float)']
     settings_default['dxmax'] = -1.0
     settings_description['dxmax'] = 'Maximum panel size'
 
@@ -86,9 +86,6 @@ def initialise(self, data, in_dict=None, restart=False):
             elif 'StaticCoupled' in data.settings:
                 aero_solver_name = data.settings['StaticCoupled']['aero_solver']
                 aero_solver_settings = data.settings['StaticCoupled']['aero_solver_settings']
-            elif 'StaticCoupled' in data.settings:
-                aero_solver_name = data.settings['StaticCoupled']['aero_solver']
-                aero_solver_settings = data.settings['StaticCoupled']['aero_solver_settings']
             elif 'DynamicCoupled' in data.settings:
                 aero_solver_name = data.settings['DynamicCoupled']['aero_solver']
                 aero_solver_settings = data.settings['DynamicCoupled']['aero_solver_settings']
@@ -121,18 +118,11 @@ def initialise(self, data, in_dict=None, restart=False):
         self.u_inf_direction = self.in_dict['u_inf_direction']
         self.dt = self.in_dict['dt']
 
-        if self.in_dict['dx1'] == -1:
-            self.dx1 = self.u_inf*self.dt
-        else:
-            self.dx1 = self.in_dict['dx1']
-
+        self.dx1 = self.in_dict['dx1']
         self.ndx1 = self.in_dict['ndx1']
         self.r = self.in_dict['r']
+        self.dxmax = self.in_dict['dxmax']
 
-        if self.in_dict['dxmax'] == -1:
-            self.dxmax = self.dx1
-        else:
-            self.dxmax = self.in_dict['dxmax']
 
     def generate(self, params):
         # Renaming for convenience
@@ -144,11 +134,12 @@ def generate(self, params):
 
         nsurf = len(zeta)
         for isurf in range(nsurf):
+            r_surf, ndx1_surf, dx1_surf, dxmax_surf = self.get_all_surface_parameters(isurf)
             M, N = zeta_star[isurf][0, :, :].shape
             for j in range(N):
                 zeta_star[isurf][:, 0, j] = zeta[isurf][:, -1, j]
                 for i in range(1, M):
-                    deltax = self.get_deltax(i, self.dx1, self.ndx1, self.r, self.dxmax)
+                    deltax = self.get_deltax(i, dx1_surf, ndx1_surf, r_surf, dxmax_surf)
                     zeta_star[isurf][:, i, j] = zeta_star[isurf][:, i - 1, j] + deltax*self.u_inf_direction
             gamma[isurf] *= 0.
             gamma_star[isurf] *= 0.
@@ -173,3 +164,23 @@ def get_deltax(i, dx1, ndx1, r, dxmax):
         deltax = min(deltax, dxmax)
 
         return deltax
+    @staticmethod
+    def get_surface_parameter(parameter, i_surf):
+        if np.isscalar(parameter):
+            return parameter
+        else:
+            return parameter[i_surf]
+
+    def get_all_surface_parameters(self,i_surf):
+        r_surf = self.get_surface_parameter(self.r, i_surf)
+        dx1_surf = self.get_surface_parameter(self.dx1, i_surf)
+        ndx1_surf = self.get_surface_parameter(self.ndx1, i_surf)
+        dxmax_surf = self.get_surface_parameter(self.dxmax, i_surf)
+        if dx1_surf == -1.:
+            dx1_surf = self.u_inf*self.dt
+        if dxmax_surf == -1.:
+            dxmax_surf = dx1_surf
+        return r_surf, ndx1_surf, dx1_surf, dxmax_surf
+
+
+
diff --git a/sharpy/linear/src/lin_aeroelastic.py b/sharpy/linear/src/lin_aeroelastic.py
index 3a1244532..e359d47c1 100644
--- a/sharpy/linear/src/lin_aeroelastic.py
+++ b/sharpy/linear/src/lin_aeroelastic.py
@@ -70,8 +70,6 @@ def __init__(self, data, custom_settings_linear=None, uvlm_block=False, chosen_t
         except KeyError:
             self.rigid_body_motions = False
 
-        print("rbm:", self.rigid_body_motions)
-
         try:
             self.use_euler = settings_here['use_euler']
         except KeyError:
@@ -104,8 +102,7 @@ def __init__(self, data, custom_settings_linear=None, uvlm_block=False, chosen_t
             rho = data.settings['StaticCoupled']['aero_solver_settings']['rho']
             if type(rho) == str:
                 rho = np.float(rho)
-            if hasattr(rho, 'value'):
-                rho = rho
+ 
         self.tsaero.rho = rho
         # --- backward compatibility
 
diff --git a/sharpy/postproc/aeroforcescalculator.py b/sharpy/postproc/aeroforcescalculator.py
index 0e44bc621..fce5468b8 100644
--- a/sharpy/postproc/aeroforcescalculator.py
+++ b/sharpy/postproc/aeroforcescalculator.py
@@ -3,7 +3,7 @@
 
 import sharpy.utils.cout_utils as cout
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.algebra as algebra
 import sharpy.aero.utils.mapping as mapping
 
@@ -54,7 +54,7 @@ class AeroForcesCalculator(BaseSolver):
     settings_default['c_ref'] = 1
     settings_description['c_ref'] = 'Reference chord'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -72,7 +72,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
         self.data = data
         self.settings = data.settings[self.solver_id]
         self.ts_max = len(self.data.structure.timestep_info)
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
         self.caller = caller
 
         self.folder = data.output_folder + '/forces/'
@@ -89,7 +89,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
 
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         if online:
             self.ts_max = len(self.data.structure.timestep_info)
diff --git a/sharpy/postproc/asymptoticstability.py b/sharpy/postproc/asymptoticstability.py
index 7fabf4c87..3565dbec7 100644
--- a/sharpy/postproc/asymptoticstability.py
+++ b/sharpy/postproc/asymptoticstability.py
@@ -2,7 +2,7 @@
 import warnings
 import numpy as np
 import scipy.linalg as sclalg
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.solver_interface import solver, BaseSolver, initialise_solver
 import sharpy.utils.cout_utils as cout
 import sharpy.utils.algebra as algebra
@@ -89,7 +89,7 @@ class AsymptoticStability(BaseSolver):
     settings_description['modes_to_plot'] = 'List of mode numbers to plot. Plots the 0, 45, 90 and 135' \
                                             'degree phases.'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -113,7 +113,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
         else:
             self.settings = custom_settings
 
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default, 
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default, 
                            options=self.settings_options, no_ctype=True)
 
         self.num_evals = self.settings['num_evals']
@@ -146,7 +146,7 @@ def run(self, **kwargs):
             eigenvectors (np.ndarray): Corresponding mode shapes
 
         """
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         # if the system is scaled, only one system can be analysed
         if self.settings['reference_velocity'] != 1. and self.data.linear.linear_system.uvlm.scaled:
diff --git a/sharpy/postproc/beamloads.py b/sharpy/postproc/beamloads.py
index b2291307a..2e113342a 100644
--- a/sharpy/postproc/beamloads.py
+++ b/sharpy/postproc/beamloads.py
@@ -2,7 +2,7 @@
 import numpy as np
 import os
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.structure.utils.xbeamlib as xbeamlib
 
 
@@ -27,7 +27,7 @@ class BeamLoads(BaseSolver):
     settings_default['output_file_name'] = 'beam_loads'
     settings_description['output_file_name'] = 'Output file name'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -44,7 +44,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
         self.caller = caller
 
         self.folder = data.output_folder + '/beam/'
@@ -53,7 +53,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
 
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         self.calculate_loads(online)
         if self.settings['csv_output']:
diff --git a/sharpy/postproc/beamplot.py b/sharpy/postproc/beamplot.py
index 5e892170e..8c63feae3 100644
--- a/sharpy/postproc/beamplot.py
+++ b/sharpy/postproc/beamplot.py
@@ -5,7 +5,7 @@
 
 import sharpy.utils.cout_utils as cout
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.algebra as algebra
 
 
@@ -45,7 +45,7 @@ class BeamPlot(BaseSolver):
     settings_default['output_rbm'] = True
     settings_description['output_rbm'] = 'Write ``csv`` file with rigid body motion data'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -64,7 +64,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
         # create folder for containing files if necessary
         self.folder = data.output_folder + '/beam/'
         if not os.path.exists(self.folder):
@@ -81,7 +81,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
 
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         self.plot(online)
         if not online:
diff --git a/sharpy/postproc/cleanup.py b/sharpy/postproc/cleanup.py
index 34312f454..019071f75 100644
--- a/sharpy/postproc/cleanup.py
+++ b/sharpy/postproc/cleanup.py
@@ -3,11 +3,8 @@
 import numpy as np
 from tvtk.api import tvtk, write_data
 
-import sharpy.utils.cout_utils as cout
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
-import sharpy.utils.algebra as algebra
-import sharpy.structure.utils.xbeamlib as xbeamlib
+import sharpy.utils.settings as settings_utils
 
 
 @solver
@@ -34,7 +31,7 @@ class Cleanup(BaseSolver):
     settings_default['remaining_steps'] = 10
     settings_description['remaining_steps'] = 'The last `remaining_steps` are not cleaned up'
 
-    table = su.SettingsTable()
+    table = settings_utils.SettingsTable()
     __doc__ += table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -48,14 +45,14 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default)
         self.caller = caller
 
     def run(self, **kwargs):
         
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         if self.settings['clean_structure']:
             self.clean(self.data.structure.timestep_info, self.settings['remaining_steps'])
diff --git a/sharpy/postproc/frequencyresponse.py b/sharpy/postproc/frequencyresponse.py
index a560a269b..6c5be1a7b 100644
--- a/sharpy/postproc/frequencyresponse.py
+++ b/sharpy/postproc/frequencyresponse.py
@@ -2,7 +2,7 @@
 import time
 import os
 import sharpy.utils.solver_interface as solver_interface
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.cout_utils as cout
 import warnings
 import sharpy.linear.src.libss as libss
@@ -92,10 +92,10 @@ class FrequencyResponse(solver_interface.BaseSolver):
     settings_default['quick_plot'] = False
     settings_description['quick_plot'] = 'Produce array of ``.png`` plots showing response. Requires matplotlib.'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
 
-    scaling_table = su.SettingsTable()
+    scaling_table = settings_utils.SettingsTable()
     __doc__ += scaling_table.generate(scaling_types, scaling_default, scaling_description,
                                       header_line='The scaling dictionary takes the following entries:')
 
@@ -119,7 +119,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
             self.settings = self.data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            self.settings_options,
@@ -133,7 +133,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
                 scaling = self.data.linear.linear_system.uvlm.sys.ScalingFacts
             else:
                 scaling = self.settings['frequency_scaling']
-                su.to_custom_types(scaling, self.scaling_types, self.scaling_default,
+                settings_utils.to_custom_types(scaling, self.scaling_types, self.scaling_default,
                                                no_ctype=True)
             self.w_to_k = scaling['length'] / scaling['speed']
         else:
@@ -166,8 +166,8 @@ def run(self, **kwargs):
               If not given, the response for the previously assembled systems and specified in ``target_system`` will
               be performed.
         """
-        online = su.set_value_or_default(kwargs, 'online', False)
-        ss = su.set_value_or_default(kwargs, 'ss', None)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
+        ss = settings_utils.set_value_or_default(kwargs, 'ss', None)
 
         if ss is None:
             ss_list = [find_target_system(self.data, system_name) for system_name in self.settings['target_system']]
diff --git a/sharpy/postproc/liftdistribution.py b/sharpy/postproc/liftdistribution.py
index 2eb03611c..d3e2c7cd0 100644
--- a/sharpy/postproc/liftdistribution.py
+++ b/sharpy/postproc/liftdistribution.py
@@ -5,9 +5,8 @@
 import sharpy.utils.cout_utils as cout
 import sharpy.utils.algebra as algebra
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.datastructures import init_matrix_structure, standalone_ctypes_pointer
-import sharpy.aero.utils.uvlmlib as uvlmlib
 import sharpy.aero.utils.mapping as mapping
 import sharpy.aero.utils.utils as aeroutils
 
@@ -38,7 +37,7 @@ class LiftDistribution(BaseSolver):
     settings_default['rho'] = 1.225
     settings_description['rho'] = 'Reference freestream density [kg/m³]'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -53,7 +52,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
         self.ts_max = len(self.data.structure.timestep_info)
         self.caller = caller
         self.folder = data.output_folder
@@ -62,7 +61,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
 
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         if not online:
             for self.ts in range(self.ts_max):
diff --git a/sharpy/postproc/pickledata.py b/sharpy/postproc/pickledata.py
index 386c0b43f..f85315828 100644
--- a/sharpy/postproc/pickledata.py
+++ b/sharpy/postproc/pickledata.py
@@ -4,10 +4,8 @@
 import h5py
 
 import sharpy
-import sharpy.utils.cout_utils as cout
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
-import sharpy.utils.h5utils as h5utils
+import sharpy.utils.settings as settings_utils
 
 
 # Define basic numerical types
@@ -35,7 +33,7 @@ class PickleData(BaseSolver):
     settings_default['stride'] = 1
     settings_description['stride'] = 'Number of steps between the execution calls when run online'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -53,7 +51,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
         else:
             self.settings = custom_settings
 
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default)
 
@@ -64,8 +62,8 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
         
-        online = su.set_value_or_default(kwargs, 'online', False)
-        solvers = su.set_value_or_default(kwargs, 'solvers', None)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
+        solvers = settings_utils.set_value_or_default(kwargs, 'solvers', None)
         
         if ((online and (self.data.ts % self.settings['stride'] == 0)) or (not online)):
             with open(self.filename, 'wb') as f:
diff --git a/sharpy/postproc/plotflowfield.py b/sharpy/postproc/plotflowfield.py
index b688e16af..f11adacf0 100644
--- a/sharpy/postproc/plotflowfield.py
+++ b/sharpy/postproc/plotflowfield.py
@@ -3,7 +3,7 @@
 from tvtk.api import tvtk, write_data
 from sharpy.utils.solver_interface import solver, BaseSolver
 import sharpy.utils.generator_interface as gen_interface
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.aero.utils.uvlmlib as uvlmlib
 import ctypes as ct
 from sharpy.utils.constants import vortex_radius_def
@@ -63,7 +63,7 @@ class PlotFlowField(BaseSolver):
     settings_default['vortex_radius'] = vortex_radius_def
     settings_description['vortex_radius'] = 'Distance below which inductions are not computed.'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
 
     def __init__(self):
@@ -78,7 +78,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            self.settings_options)
@@ -184,7 +184,7 @@ def output_velocity_field(self, ts):
 
     def run(self, **kwargs):
 
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         if online:
             if divmod(self.data.ts, self.settings['stride'])[1] == 0:
diff --git a/sharpy/postproc/savedata.py b/sharpy/postproc/savedata.py
index 8d06fa85a..50719e0bf 100644
--- a/sharpy/postproc/savedata.py
+++ b/sharpy/postproc/savedata.py
@@ -4,7 +4,7 @@
 import sharpy
 import sharpy.utils.cout_utils as cout
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.h5utils as h5utils
 from sharpy.presharpy.presharpy import PreSharpy
 
@@ -94,7 +94,7 @@ class SaveData(BaseSolver):
     settings_default['stride'] = 1
     settings_description['stride'] = 'Number of steps between the execution calls when run online'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description,
                                        settings_options=settings_options)
 
@@ -120,7 +120,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
         else:
             self.settings = custom_settings
 
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            options=self.settings_options)
@@ -197,7 +197,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
 
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         # Use the following statement in case the ct types are not defined and
         # you need them on uvlm3d
diff --git a/sharpy/postproc/saveparametriccase.py b/sharpy/postproc/saveparametriccase.py
index a61508a47..89db7bf61 100644
--- a/sharpy/postproc/saveparametriccase.py
+++ b/sharpy/postproc/saveparametriccase.py
@@ -1,5 +1,5 @@
 from sharpy.utils.solver_interface import solver, BaseSolver, initialise_solver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import configobj
 import os
 import sharpy.utils.cout_utils as cout
@@ -71,7 +71,7 @@ class SaveParametricCase(BaseSolver):
                                                    'the M, C, K matrices. The setting ``save_pmor_items`` ' \
                                                    'should be set to `on`'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -88,7 +88,7 @@ def initialise(self, data, custom_settings=None, restart=False):
         else:
             self.settings = custom_settings
 
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default)
 
@@ -97,8 +97,8 @@ def initialise(self, data, custom_settings=None, restart=False):
 
     def run(self, **kwargs):
         
-        online = su.set_value_or_default(kwargs, 'online', False)
-        restart = su.set_value_or_default(kwargs, 'restart', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
+        restart = settings_utils.set_value_or_default(kwargs, 'restart', False)
 
         config = configobj.ConfigObj()
         file_name = self.folder + '/' + self.data.settings['SHARPy']['case'] + '.pmor.sharpy'
diff --git a/sharpy/postproc/stabilityderivatives.py b/sharpy/postproc/stabilityderivatives.py
index 8da5857c9..30d5d43c0 100644
--- a/sharpy/postproc/stabilityderivatives.py
+++ b/sharpy/postproc/stabilityderivatives.py
@@ -2,9 +2,8 @@
 import os
 import numpy as np
 import sharpy.utils.cout_utils as cout
-import sharpy.utils.algebra as algebra
-import sharpy.utils.settings as su
-from sharpy.linear.utils.derivatives import Derivatives, DerivativeSet
+import sharpy.utils.settings as settings_utils
+from sharpy.linear.utils.derivatives import Derivatives
 
 
 @solver_interface.solver
@@ -41,7 +40,7 @@ class StabilityDerivatives(solver_interface.BaseSolver):
     settings_default['c_ref'] = 1.
     settings_description['c_ref'] = 'Reference chord'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -66,7 +65,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
         else:
             self.settings = self.data.settings[self.solver_id]
 
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            options=self.settings_options,
@@ -107,7 +106,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
         
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         # TODO: consider running all required solvers inside this one to keep the correct settings
         # i.e: run Modal, Linear Assembly
diff --git a/sharpy/postproc/stallcheck.py b/sharpy/postproc/stallcheck.py
index f88fc29fa..e75fbff29 100644
--- a/sharpy/postproc/stallcheck.py
+++ b/sharpy/postproc/stallcheck.py
@@ -1,7 +1,7 @@
 import numpy as np
 import sharpy.utils.cout_utils as cout
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.datastructures import init_matrix_structure, standalone_ctypes_pointer
 import sharpy.aero.utils.uvlmlib as uvlmlib
 
@@ -43,7 +43,7 @@ class StallCheck(BaseSolver):
     settings_default['output_degrees'] = False
     settings_description['output_degrees'] = 'Output incidence angles in degrees vs radians'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -61,13 +61,13 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
         self.ts_max = len(self.data.structure.timestep_info)
         self.caller = caller
 
     def run(self, **kwargs):
     
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         if not online:
             for self.ts in range(self.ts_max):
diff --git a/sharpy/postproc/udpout.py b/sharpy/postproc/udpout.py
index fb5d52639..6751264a4 100644
--- a/sharpy/postproc/udpout.py
+++ b/sharpy/postproc/udpout.py
@@ -1,5 +1,5 @@
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.io.network_interface as network_interface
 import sharpy.utils.cout_utils as cout
 
@@ -33,7 +33,7 @@ class UDPout(BaseSolver):
     del settings_default['send_output_to_all_clients']
     del settings_description['send_output_to_all_clients']
 
-    table = su.SettingsTable()
+    table = settings_utils.SettingsTable()
     __doc__ += table.generate(settings_types, settings_default, settings_description,
                               header_line='This post-processor takes in the following settings, for a more '
                                           'detailed description see '
@@ -71,7 +71,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
         
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         if online:
             self.set_of_variables.get_value(self.data)
diff --git a/sharpy/postproc/writevariablestime.py b/sharpy/postproc/writevariablestime.py
index 5be78a963..fe532d8df 100644
--- a/sharpy/postproc/writevariablestime.py
+++ b/sharpy/postproc/writevariablestime.py
@@ -1,7 +1,7 @@
 import os
 import numpy as np
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 
 
 @solver
@@ -92,7 +92,7 @@ class WriteVariablesTime(BaseSolver):
     settings_default['vel_field_points'] = np.array([0., 0., 0.])
     settings_description['vel_field_points'] = 'List of coordinates of the control points as x1, y1, z1, x2, y2, z2 ...'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -111,7 +111,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         self.folder = data.output_folder + '/WriteVariablesTime/'
         if not os.path.isdir(self.folder):
@@ -195,7 +195,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
     
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = settings_utils.set_value_or_default(kwargs, 'online', False)
 
         if online:
             self.data = self.write(-1)
diff --git a/sharpy/solvers/aerogridloader.py b/sharpy/solvers/aerogridloader.py
index 99fe104d3..bc8f9a8df 100644
--- a/sharpy/solvers/aerogridloader.py
+++ b/sharpy/solvers/aerogridloader.py
@@ -3,7 +3,7 @@
 
 from sharpy.utils.solver_interface import solver, BaseSolver
 import sharpy.aero.models.aerogrid as aerogrid
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.h5utils as h5utils
 import sharpy.utils.generator_interface as gen_interface
 
@@ -62,7 +62,7 @@ class AerogridLoader(BaseSolver):
     settings_default['freestream_dir'] = [1.0, 0.0, 0.0]
     settings_description['freestream_dir'] = 'Free stream flow direction'
 
-    settings_types['mstar'] = 'int'
+    settings_types['mstar'] = ['int', 'list(int)']
     settings_default['mstar'] = 10
     settings_description['mstar'] = 'Number of chordwise wake panels'
 
@@ -84,7 +84,7 @@ class AerogridLoader(BaseSolver):
     settings_default['wake_shape_generator_input'] = dict()
     settings_description['wake_shape_generator_input'] = 'Dictionary of inputs needed by the wake shape generator'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description,
                                        settings_options=settings_options)
 
@@ -105,7 +105,7 @@ def initialise(self, data, restart=False):
         self.settings = data.settings[self.solver_id]
 
         # init settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default, options=self.settings_options)
 
diff --git a/sharpy/solvers/beamloader.py b/sharpy/solvers/beamloader.py
index ecdb2f3e2..a30ee6086 100644
--- a/sharpy/solvers/beamloader.py
+++ b/sharpy/solvers/beamloader.py
@@ -1,9 +1,8 @@
 import h5py as h5
-import numpy as np
 
 from sharpy.utils.solver_interface import solver, BaseSolver
 import sharpy.structure.models.beam as beam
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.h5utils as h5utils
 import os
 
@@ -60,7 +59,7 @@ class BeamLoader(BaseSolver):
     settings_default['for_pos'] = [0., 0, 0]
     settings_description['for_pos'] = 'Initial position of the A FoR.'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -81,7 +80,7 @@ def initialise(self, data, restart=False):
         self.settings = data.settings[self.solver_id]
 
         # init settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         # read input files (fem and dyn)
         self.read_files()
diff --git a/sharpy/solvers/dynamiccoupled.py b/sharpy/solvers/dynamiccoupled.py
index 1bb97c3ac..72d107467 100644
--- a/sharpy/solvers/dynamiccoupled.py
+++ b/sharpy/solvers/dynamiccoupled.py
@@ -13,7 +13,7 @@
 import sharpy.utils.solver_interface as solver_interface
 import sharpy.utils.controller_interface as controller_interface
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.algebra as algebra
 import sharpy.utils.exceptions as exc
 import sharpy.io.network_interface as network_interface
@@ -178,7 +178,7 @@ class DynamicCoupled(BaseSolver):
                                                  'The dictionary values are dictionaries with the settings ' \
                                                  'needed by each generator.'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
 
     def __init__(self):
@@ -254,7 +254,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            options=self.settings_options)
@@ -341,7 +341,8 @@ def initialise(self, data, custom_settings=None, restart=False):
                                                      aero=self.data.aero,
                                                      structure=self.data.structure,
                                                      rho=self.settings['aero_solver_settings']['rho'],
-                                                     vortex_radius=self.settings['aero_solver_settings']['vortex_radius'])
+                                                     vortex_radius=self.settings['aero_solver_settings']['vortex_radius'],
+                                                     output_folder = self.data.output_folder)
 
         # check for empty dictionary
         if self.settings['network_settings']:
@@ -442,7 +443,7 @@ def run(self, **kwargs):
         Run the time stepping procedure with controllers and postprocessors
         included.
         """
-        solvers = su.set_value_or_default(kwargs, 'solvers', None)
+        solvers = settings_utils.set_value_or_default(kwargs, 'solvers', None)
         if self.network_loader is not None:
             self.set_of_variables = self.network_loader.get_inout_variables()
 
@@ -817,7 +818,8 @@ def map_forces(self, aero_kstep, structural_kstep, unsteady_forces_coeff=1.0):
             struct_forces = \
                 self.correct_forces_generator.generate(aero_kstep=aero_kstep,
                                                        structural_kstep=structural_kstep,
-                                                       struct_forces=struct_forces)
+                                                       struct_forces=struct_forces,
+                                                       ts=self.data.ts)
 
         aero_kstep.aero_steady_forces_beam_dof = struct_forces
         structural_kstep.postproc_node['aero_steady_forces'] = struct_forces
diff --git a/sharpy/solvers/dynamicuvlm.py b/sharpy/solvers/dynamicuvlm.py
index c0a2be909..a813677f3 100644
--- a/sharpy/solvers/dynamicuvlm.py
+++ b/sharpy/solvers/dynamicuvlm.py
@@ -3,13 +3,9 @@
 
 N Goizueta Jan 19
 """
-
-import ctypes as ct
-import numpy as np
-
 import sharpy.utils.solver_interface as solver_interface
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.cout_utils as cout
 
 @solver
@@ -74,7 +70,7 @@ class DynamicUVLM(BaseSolver):
     settings_default['postprocessors_settings'] = dict()
     settings_description['postprocessors_settings'] = 'Dictionary with the applicable settings for every ``psotprocessor``. Every ``postprocessor`` needs its entry, even if empty'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
 
@@ -98,7 +94,7 @@ def initialise(self, data, custom_settings=None, restart=False):
         else:
             self.settings = custom_settings
 
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
         self.dt = self.settings['dt']
         self.print_info = self.settings['print_info']
 
diff --git a/sharpy/solvers/initialaeroelasticloader.py b/sharpy/solvers/initialaeroelasticloader.py
index 41064711b..5adabf534 100644
--- a/sharpy/solvers/initialaeroelasticloader.py
+++ b/sharpy/solvers/initialaeroelasticloader.py
@@ -1,9 +1,5 @@
-import numpy as np
-import h5py as h5
-
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
-import sharpy.utils.algebra as algebra
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.h5utils as h5utils
 import sharpy.utils.exceptions as exceptions
 
@@ -33,7 +29,7 @@ class InitialAeroelasticLoader(BaseSolver):
     settings_default['generate_aero'] = False
     settings_description['generate_aero'] = 'Generate the aerodynamics grids from scratch'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -48,7 +44,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            no_ctype=True)
@@ -57,8 +53,8 @@ def initialise(self, data, custom_settings=None, restart=False):
 
     def run(self, **kwargs):
 
-        aero_step = su.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
-        structural_step = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        aero_step = settings_utils.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
+        structural_step = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
         
         # Copy structural information
         attributes = ['pos', 'pos_dot', 'pos_ddot',
diff --git a/sharpy/solvers/lindynamicsim.py b/sharpy/solvers/lindynamicsim.py
index 27e5a5a8f..8538e0bab 100644
--- a/sharpy/solvers/lindynamicsim.py
+++ b/sharpy/solvers/lindynamicsim.py
@@ -2,7 +2,7 @@
 import os
 import h5py as h5
 from sharpy.utils.solver_interface import solver, BaseSolver, initialise_solver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.linear.src.libss as libss
 import scipy.linalg as sclalg
 import sharpy.utils.h5utils as h5utils
@@ -73,7 +73,7 @@ class LinDynamicSim(BaseSolver):
     settings_types['postprocessors_settings'] = 'dict'
     settings_default['postprocessors_settings'] = dict()
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -95,7 +95,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = custom_settings
         else:
             self.settings = data.settings[self.solver_id]
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default, no_ctype=True)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default, no_ctype=True)
 
         # Read initial state and input data and store in dictionary
         self.read_files()
diff --git a/sharpy/solvers/linearassembler.py b/sharpy/solvers/linearassembler.py
index b9dce02bd..69f4908f4 100644
--- a/sharpy/solvers/linearassembler.py
+++ b/sharpy/solvers/linearassembler.py
@@ -5,8 +5,7 @@
 from sharpy.utils.solver_interface import solver, BaseSolver
 
 import sharpy.linear.utils.ss_interface as ss_interface
-import sharpy.utils.settings as su
-import sharpy.utils.h5utils as h5
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.cout_utils as cout
 
 
@@ -105,7 +104,7 @@ class LinearAssembler(BaseSolver):
     settings_default['recover_accelerations'] = False
     settings_description['recover_accelerations'] = 'Recover structural system accelerations as additional outputs.'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
 
     def __init__(self):
@@ -123,7 +122,7 @@ def initialise(self, data, custom_settings=None, restart=False):
 
         else:
             self.settings = data.settings[self.solver_id]
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            options=self.settings_options,
diff --git a/sharpy/solvers/modal.py b/sharpy/solvers/modal.py
index 76a462038..3fc053a8e 100644
--- a/sharpy/solvers/modal.py
+++ b/sharpy/solvers/modal.py
@@ -6,7 +6,7 @@
 import warnings
 import sharpy.structure.utils.xbeamlib as xbeamlib
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.algebra as algebra
 import sharpy.utils.cout_utils as cout
 import sharpy.structure.utils.modalutils as modalutils
@@ -95,7 +95,7 @@ class Modal(BaseSolver):
     settings_default['rigid_modes_cg'] = False
     settings_description['rigid_modes_cg'] = 'Not implemente yet'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -115,7 +115,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default)
 
@@ -236,7 +236,7 @@ def run(self, **kwargs):
                 # full_matrix_settings = self.data.settings['DynamicCoupled']['structural_solver_settings']
             import sharpy.solvers._basestructural as basestructuralsolver
             full_matrix_settings = basestructuralsolver._BaseStructural().settings_default
-            su.to_custom_types(full_matrix_settings, basestructuralsolver._BaseStructural().settings_types, full_matrix_settings)
+            settings_utils.to_custom_types(full_matrix_settings, basestructuralsolver._BaseStructural().settings_types, full_matrix_settings)
 
 
             # Obtain the tangent mass, damping and stiffness matrices
diff --git a/sharpy/solvers/noaero.py b/sharpy/solvers/noaero.py
index ae1760b17..5a18a68d9 100644
--- a/sharpy/solvers/noaero.py
+++ b/sharpy/solvers/noaero.py
@@ -1,4 +1,4 @@
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.solver_interface import solver, BaseSolver
 
 
@@ -28,7 +28,7 @@ class NoAero(BaseSolver):
     settings_default['update_grid'] = True
     settings_description['update_grid'] = 'Update aerodynamic grid as the structure deforms.'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -41,7 +41,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            no_ctype=True)
@@ -51,12 +51,12 @@ def initialise(self, data, custom_settings=None, restart=False):
 
     def run(self, **kwargs):
 
-        aero_tstep = su.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
-        structure_tstep = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
-        convect_wake = su.set_value_or_default(kwargs, 'convect_wake', True)
-        dt= su.set_value_or_default(kwargs, 'dt', self.settings['dt'])                                                                                                                             
-        t = su.set_value_or_default(kwargs, 't', self.data.ts*dt)
-        unsteady_contribution = su.set_value_or_default(kwargs, 'unsteady_contribution', False)
+        aero_tstep = settings_utils.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
+        structure_tstep = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        convect_wake = settings_utils.set_value_or_default(kwargs, 'convect_wake', True)
+        dt= settings_utils.set_value_or_default(kwargs, 'dt', self.settings['dt'])                                                                                                                             
+        t = settings_utils.set_value_or_default(kwargs, 't', self.data.ts*dt)
+        unsteady_contribution = settings_utils.set_value_or_default(kwargs, 'unsteady_contribution', False)
 
         # generate the wake because the solid shape might change
         self.data.aero.wake_shape_generator.generate({'zeta': aero_tstep.zeta,
diff --git a/sharpy/solvers/nonlineardynamic.py b/sharpy/solvers/nonlineardynamic.py
index 97200336e..953305ffe 100644
--- a/sharpy/solvers/nonlineardynamic.py
+++ b/sharpy/solvers/nonlineardynamic.py
@@ -2,13 +2,10 @@
 @modified   Alfonso del Carre
 """
 
-import ctypes as ct
-import numpy as np
-
 import sharpy.structure.utils.xbeamlib as xbeamlib
 from sharpy.utils.settings import str2bool
-from sharpy.utils.solver_interface import solver, BaseSolver, solver_from_string
-import sharpy.utils.settings as su
+from sharpy.utils.solver_interface import solver, solver_from_string
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.cout_utils as cout
 
 _BaseStructural = solver_from_string('_BaseStructural')
@@ -38,7 +35,7 @@ class NonLinearDynamic(_BaseStructural):
     settings_types['gravity_dir'] = 'list(float)'
     settings_default['gravity_dir'] = [0, 0, 1]
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -48,7 +45,7 @@ def __init__(self):
     def initialise(self, data, restart=False):
         self.data = data
         self.settings = data.settings[self.solver_id]
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         # load info from dyn dictionary
         self.data.structure.add_unsteady_information(self.data.structure.dyn_dict, self.settings['num_steps'])
diff --git a/sharpy/solvers/nonlineardynamiccoupledstep.py b/sharpy/solvers/nonlineardynamiccoupledstep.py
index 80da34353..621686ba5 100644
--- a/sharpy/solvers/nonlineardynamiccoupledstep.py
+++ b/sharpy/solvers/nonlineardynamiccoupledstep.py
@@ -2,15 +2,12 @@
 @modified   Alfonso del Carre
 """
 
-import ctypes as ct
 import numpy as np
 
 import sharpy.structure.utils.xbeamlib as xbeamlib
 from sharpy.utils.settings import str2bool
 from sharpy.utils.solver_interface import solver, BaseSolver, solver_from_string
-import sharpy.utils.settings as su
-import sharpy.utils.algebra as algebra
-import sharpy.utils.cout_utils as cout
+import sharpy.utils.settings as settings_utils
 
 
 _BaseStructural = solver_from_string('_BaseStructural')
@@ -50,7 +47,7 @@ class NonLinearDynamicCoupledStep(_BaseStructural):
     settings_default['relaxation_factor'] = 0.3
     settings_description['relaxation factor'] = 'Relaxation factor'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -63,7 +60,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         # load info from dyn dictionary
         self.data.structure.add_unsteady_information(self.data.structure.dyn_dict, self.settings['num_steps'])
@@ -79,10 +76,10 @@ def initialise(self, data, custom_settings=None, restart=False):
 
     def run(self, **kwargs):
 
-        structural_step = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        structural_step = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
         # TODO: previous_structural_step never used
-        previous_structural_step = su.set_value_or_default(kwargs, 'previous_structural_step', self.data.structure.timestep_info[-1])
-        dt= su.set_value_or_default(kwargs, 'dt', self.settings['dt'])
+        previous_structural_step = settings_utils.set_value_or_default(kwargs, 'previous_structural_step', self.data.structure.timestep_info[-1])
+        dt= settings_utils.set_value_or_default(kwargs, 'dt', self.settings['dt'])
 
         xbeamlib.xbeam_step_couplednlndyn(self.data.structure,
                                           self.settings,
diff --git a/sharpy/solvers/nonlineardynamicmultibody.py b/sharpy/solvers/nonlineardynamicmultibody.py
index fb45e24ac..bcbd4be7d 100644
--- a/sharpy/solvers/nonlineardynamicmultibody.py
+++ b/sharpy/solvers/nonlineardynamicmultibody.py
@@ -3,7 +3,7 @@
 import os
 
 from sharpy.utils.solver_interface import solver, BaseSolver, solver_from_string
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.solver_interface as solver_interface
 
 import sharpy.utils.cout_utils as cout
@@ -62,7 +62,7 @@ class NonLinearDynamicMultibody(_BaseStructural):
     settings_default['zero_ini_dot_ddot'] = False
     settings_description['zero_ini_dot_ddot'] = 'Set to zero the position and crv derivatives at the first time step'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -93,7 +93,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            self.settings_options,
@@ -350,8 +350,8 @@ def write_lm_cond_num(self, iteration, Lambda, Lambda_dot, Lambda_ddot, cond_num
 
 
     def run(self, **kwargs):
-        structural_step = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
-        dt= su.set_value_or_default(kwargs, 'dt', self.settings['dt'])
+        structural_step = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        dt= settings_utils.set_value_or_default(kwargs, 'dt', self.settings['dt'])
 
         if structural_step.mb_dict is not None:
             MBdict = structural_step.mb_dict
@@ -399,7 +399,6 @@ def run(self, **kwargs):
         old_Dq = 1.0
         LM_old_Dq = 1.0
 
-        converged = False
         skip_step = False
         for iteration in range(self.settings['max_iterations']):
             # Check if the maximum of iterations has been reached
@@ -467,7 +466,6 @@ def run(self, **kwargs):
             # Evaluate convergence
             res = np.max(np.abs(Dq[0:self.sys_size]))/old_Dq
             if np.isnan(res):
-                print(old_Dq)
                 if self.settings['allow_skip_step']:
                     skip_step = True
                     cout.cout_wrap("Skipping step", 3)
@@ -480,8 +478,6 @@ def run(self, **kwargs):
                 LM_res = 0.0
 
             if (res < self.settings['min_delta']) and (LM_res < self.settings['min_delta']):
-            # if (res < self.settings['min_delta']):
-                converged = True
                 break
 
         Lambda, Lambda_dot = mb.state2disp_and_accel(q, dqdt, dqddt, MB_beam, MB_tstep, num_LM_eq)
diff --git a/sharpy/solvers/nonlineardynamicprescribedstep.py b/sharpy/solvers/nonlineardynamicprescribedstep.py
index f68c3762d..535a5d225 100644
--- a/sharpy/solvers/nonlineardynamicprescribedstep.py
+++ b/sharpy/solvers/nonlineardynamicprescribedstep.py
@@ -2,15 +2,9 @@
 @modified   Alfonso del Carre
 """
 
-import ctypes as ct
-import numpy as np
-
 import sharpy.structure.utils.xbeamlib as xbeamlib
-from sharpy.utils.settings import str2bool
 from sharpy.utils.solver_interface import solver, BaseSolver, solver_from_string
-import sharpy.utils.settings as su
-import sharpy.utils.cout_utils as cout
-import sharpy.utils.algebra as algebra
+import sharpy.utils.settings as settings_utils
 
 _BaseStructural = solver_from_string('_BaseStructural')
 
@@ -34,7 +28,7 @@ class NonLinearDynamicPrescribedStep(_BaseStructural):
     settings_default = _BaseStructural.settings_default.copy()
     settings_description = _BaseStructural.settings_description.copy()
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -47,15 +41,15 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         # load info from dyn dictionary
         self.data.structure.add_unsteady_information(self.data.structure.dyn_dict, self.settings['num_steps'])
 
 
     def run(self, **kwargs):
-        structural_step = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
-        dt = su.set_value_or_default(kwargs, 'dt', self.settings['dt'])  
+        structural_step = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        dt = settings_utils.set_value_or_default(kwargs, 'dt', self.settings['dt'])  
 
         if self.data.ts > 0:
             try:
diff --git a/sharpy/solvers/nonlinearstatic.py b/sharpy/solvers/nonlinearstatic.py
index 0575391b2..784608c61 100644
--- a/sharpy/solvers/nonlinearstatic.py
+++ b/sharpy/solvers/nonlinearstatic.py
@@ -1,10 +1,8 @@
 import numpy as np
 
 import sharpy.structure.utils.xbeamlib as xbeamlib
-import sharpy.utils.cout_utils as cout
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.solver_interface import solver, BaseSolver, solver_from_string
-import sharpy.utils.algebra as algebra
 
 _BaseStructural = solver_from_string('_BaseStructural')
 
@@ -34,7 +32,7 @@ class NonLinearStatic(_BaseStructural):
     settings_types['initial_velocity'] = 'list(float)'
     settings_default['initial_velocity'] = np.array([0., 0., 0., 0., 0., 0.])
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -47,7 +45,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default, no_ctype=True)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default, no_ctype=True)
 
     def run(self, **kwargs):
         self.data.structure.timestep_info[self.data.ts].for_pos[0:3] = self.settings['initial_position']
diff --git a/sharpy/solvers/nostructural.py b/sharpy/solvers/nostructural.py
index 661b9fdbf..6d5b11d0f 100644
--- a/sharpy/solvers/nostructural.py
+++ b/sharpy/solvers/nostructural.py
@@ -1,7 +1,6 @@
 import numpy as np
 
-import sharpy.utils.cout_utils as cout
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.solver_interface import solver, BaseSolver, solver_from_string
 
 _BaseStructural = solver_from_string('_BaseStructural')
@@ -29,7 +28,7 @@ class NoStructural(_BaseStructural):
     settings_types['initial_position'] = 'list(float)'
     settings_default['initial_position'] = np.array([0.0, 0.0, 0.0])
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -42,7 +41,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
     def run(self, **kwargs):
         self.data.structure.timestep_info[self.data.ts].for_pos[0:3] = self.settings['initial_position']
diff --git a/sharpy/solvers/prescribeduvlm.py b/sharpy/solvers/prescribeduvlm.py
index c02340ef8..cbd7acdb8 100644
--- a/sharpy/solvers/prescribeduvlm.py
+++ b/sharpy/solvers/prescribeduvlm.py
@@ -1,8 +1,4 @@
-import ctypes as ct
-import numpy as np
-
-import sharpy.utils.algebra as algebra
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.solver_interface import solver, BaseSolver
 import sharpy.utils.solver_interface as solver_interface
 import sharpy.utils.cout_utils as cout
@@ -70,7 +66,7 @@ class PrescribedUvlm(BaseSolver):
     settings_default['vortex_radius_wake_ind'] = vortex_radius_def
     settings_description['vortex_radius_wake_ind'] = 'Distance between points below which induction is not computed in the wake convection'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -88,7 +84,7 @@ def __init__(self):
     def initialise(self, data, restart=False):
         self.data = data
         self.settings = data.settings[self.solver_id]
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
         self.dt = self.settings['dt']
 
         self.aero_solver = solver_interface.initialise_solver(self.settings['aero_solver'])
diff --git a/sharpy/solvers/rigiddynamiccoupledstep.py b/sharpy/solvers/rigiddynamiccoupledstep.py
index f39c90eb9..74393a437 100644
--- a/sharpy/solvers/rigiddynamiccoupledstep.py
+++ b/sharpy/solvers/rigiddynamiccoupledstep.py
@@ -1,12 +1,6 @@
-import ctypes as ct
-import numpy as np
-
 import sharpy.structure.utils.xbeamlib as xbeamlib
-from sharpy.utils.settings import str2bool
 from sharpy.utils.solver_interface import solver, BaseSolver, solver_from_string
-import sharpy.utils.settings as su
-import sharpy.utils.algebra as algebra
-import sharpy.utils.cout_utils as cout
+import sharpy.utils.settings as settings_utils
 
 
 _BaseStructural = solver_from_string('_BaseStructural')
@@ -37,7 +31,7 @@ class RigidDynamicCoupledStep(_BaseStructural):
     settings_default['relaxation_factor'] = 0.3
     settings_description['relaxation factor'] = 'Relaxation factor'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -50,7 +44,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         # load info from dyn dictionary
         self.data.structure.add_unsteady_information(self.data.structure.dyn_dict, self.settings['num_steps'])
@@ -60,10 +54,10 @@ def initialise(self, data, custom_settings=None, restart=False):
 
     def run(self, **kwargs):
 
-        structural_step = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        structural_step = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
         # TODO: previous_structural_step never used
-        previous_structural_step = su.set_value_or_default(kwargs, 'previous_structural_step', self.data.structure.timestep_info[-1])
-        dt= su.set_value_or_default(kwargs, 'dt', self.settings['dt'])  
+        previous_structural_step = settings_utils.set_value_or_default(kwargs, 'previous_structural_step', self.data.structure.timestep_info[-1])
+        dt= settings_utils.set_value_or_default(kwargs, 'dt', self.settings['dt'])  
 
         xbeamlib.xbeam_step_coupledrigid(self.data.structure,
                                          self.settings,
diff --git a/sharpy/solvers/rigiddynamicprescribedstep.py b/sharpy/solvers/rigiddynamicprescribedstep.py
index f13764b57..ddff02f6d 100644
--- a/sharpy/solvers/rigiddynamicprescribedstep.py
+++ b/sharpy/solvers/rigiddynamicprescribedstep.py
@@ -2,16 +2,11 @@
 @modified   Alfonso del Carre
 """
 
-import ctypes as ct
 import numpy as np
 
 import sharpy.structure.utils.xbeamlib as xbeamlib
-# from sharpy.utils.settings import str2bool
-from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
-import sharpy.utils.cout_utils as cout
-import sharpy.utils.algebra as algebra
-from sharpy.utils.solver_interface import solver_from_string
+from sharpy.utils.solver_interface import solver, BaseSolver, solver_from_string
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.algebra as algebra
 
 _BaseStructural = solver_from_string('_BaseStructural')
@@ -33,7 +28,7 @@ class RigidDynamicPrescribedStep(BaseSolver):
     settings_default['num_steps'] = 500
     settings_description['num_steps'] = 'Number of timesteps to be run'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -47,14 +42,14 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         # load info from dyn dictionary
         self.data.structure.add_unsteady_information(self.data.structure.dyn_dict, self.settings['num_steps'])
 
     def run(self, **kwargs):
-        structural_step = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
-        dt= su.set_value_or_default(kwargs, 'dt', self.settings['dt'])  
+        structural_step = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        dt= settings_utils.set_value_or_default(kwargs, 'dt', self.settings['dt'])  
 
         if self.data.ts > 0:
             try:
diff --git a/sharpy/solvers/staticcoupled.py b/sharpy/solvers/staticcoupled.py
index 415293ec7..8fcb8a343 100644
--- a/sharpy/solvers/staticcoupled.py
+++ b/sharpy/solvers/staticcoupled.py
@@ -1,12 +1,11 @@
-import ctypes as ct
 import sys
 import numpy as np
 
 import sharpy.aero.utils.mapping as mapping
 import sharpy.utils.cout_utils as cout
-import sharpy.utils.solver_interface as solver_interface
-from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+from sharpy.utils.solver_interface import solver, BaseSolver, initialise_solver
+
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.algebra as algebra
 import sharpy.utils.generator_interface as gen_interface
 
@@ -76,7 +75,7 @@ class StaticCoupled(BaseSolver):
                                                  'The dictionary values are dictionaries with the settings ' \
                                                  'needed by each generator.'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
 
     def __init__(self):
@@ -102,7 +101,7 @@ def initialise(self, data, input_dict=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = input_dict
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            options=self.settings_options,
@@ -110,9 +109,9 @@ def initialise(self, data, input_dict=None, restart=False):
 
         self.print_info = self.settings['print_info']
 
-        self.structural_solver = solver_interface.initialise_solver(self.settings['structural_solver'])
+        self.structural_solver = initialise_solver(self.settings['structural_solver'])
         self.structural_solver.initialise(self.data, self.settings['structural_solver_settings'], restart=restart)
-        self.aero_solver = solver_interface.initialise_solver(self.settings['aero_solver'])
+        self.aero_solver = initialise_solver(self.settings['aero_solver'])
         self.aero_solver.initialise(self.structural_solver.data, self.settings['aero_solver_settings'], restart=restart)
         self.data = self.aero_solver.data
 
@@ -131,7 +130,8 @@ def initialise(self, data, input_dict=None, restart=False):
                                                      aero=self.data.aero,
                                                      structure=self.data.structure,
                                                      rho=self.settings['aero_solver_settings']['rho'],
-                                                     vortex_radius=self.settings['aero_solver_settings']['vortex_radius'])
+                                                     vortex_radius=self.settings['aero_solver_settings']['vortex_radius'],
+                                                     output_folder = self.data.output_folder)
 
         # initialise runtime generators
         self.runtime_generators = dict()
@@ -195,7 +195,8 @@ def run(self, **kwargs):
                     struct_forces = \
                         self.correct_forces_generator.generate(aero_kstep=self.data.aero.timestep_info[self.data.ts],
                                                                structural_kstep=self.data.structure.timestep_info[self.data.ts],
-                                                               struct_forces=struct_forces)
+                                                               struct_forces=struct_forces,
+                                                               ts=0)
                 self.data.aero.timestep_info[self.data.ts].aero_steady_forces_beam_dof = struct_forces
                 self.data.structure.timestep_info[self.data.ts].postproc_node['aero_steady_forces'] = struct_forces  # B
                 
diff --git a/sharpy/solvers/statictrim.py b/sharpy/solvers/statictrim.py
index d98881d0a..a76947edf 100644
--- a/sharpy/solvers/statictrim.py
+++ b/sharpy/solvers/statictrim.py
@@ -3,7 +3,7 @@
 import sharpy.utils.cout_utils as cout
 import sharpy.utils.solver_interface as solver_interface
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import os
 
 
@@ -54,7 +54,7 @@ class StaticTrim(BaseSolver):
     settings_default['m_tolerance'] = 0.01
     settings_description['m_tolerance'] = 'Tolerance in pitching moment'
 
-    settings_types['tail_cs_index'] = 'int'
+    settings_types['tail_cs_index'] = ['int', 'list(int)']
     settings_default['tail_cs_index'] = 0
     settings_description['tail_cs_index'] = 'Index of control surfaces that move to achieve trim'
 
@@ -90,7 +90,7 @@ class StaticTrim(BaseSolver):
     settings_default['save_info'] = False
     settings_description['save_info'] = 'Save trim results to text file'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -117,7 +117,7 @@ def __init__(self):
     def initialise(self, data, restart=False):
         self.data = data
         self.settings = data.settings[self.solver_id]
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         self.solver = solver_interface.initialise_solver(self.settings['solver'])
         self.solver.initialise(self.data, self.settings['solver_settings'], restart=restart)
diff --git a/sharpy/solvers/staticuvlm.py b/sharpy/solvers/staticuvlm.py
index 61eb21e0b..66fcd7d35 100644
--- a/sharpy/solvers/staticuvlm.py
+++ b/sharpy/solvers/staticuvlm.py
@@ -1,10 +1,6 @@
-import ctypes as ct
-import numpy as np
 
-import sharpy.utils.algebra as algebra
 import sharpy.aero.utils.uvlmlib as uvlmlib
-import sharpy.utils.cout_utils as cout
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.solver_interface import solver, BaseSolver
 import sharpy.utils.generator_interface as gen_interface
 from sharpy.utils.constants import vortex_radius_def
@@ -115,7 +111,7 @@ class StaticUvlm(BaseSolver):
     settings_default['map_forces_on_struct'] = False
     settings_description['map_forces_on_struct'] = 'Maps the forces on the structure at the end of the timestep. Only usefull if the solver is used outside StaticCoupled'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -130,7 +126,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default, no_ctype=True)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default, no_ctype=True)
 
         self.update_step()
 
@@ -158,10 +154,10 @@ def update_custom_grid(self, structure_tstep, aero_tstep):
 
     def run(self, **kwargs):
 
-        aero_tstep = su.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
-        structure_tstep = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
-        dt = su.set_value_or_default(kwargs, 'dt', self.settings['rollup_dt'])
-        t = su.set_value_or_default(kwargs, 't', self.data.ts*dt)
+        aero_tstep = settings_utils.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
+        structure_tstep = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        dt = settings_utils.set_value_or_default(kwargs, 'dt', self.settings['rollup_dt'])
+        t = settings_utils.set_value_or_default(kwargs, 't', self.data.ts*dt)
 
         unsteady_contribution = False
         convect_wake = False
diff --git a/sharpy/solvers/steplinearuvlm.py b/sharpy/solvers/steplinearuvlm.py
index af2edaa3c..047773b4f 100644
--- a/sharpy/solvers/steplinearuvlm.py
+++ b/sharpy/solvers/steplinearuvlm.py
@@ -5,7 +5,7 @@
 """
 from sharpy.utils.solver_interface import BaseSolver, solver
 import numpy as np
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.generator_interface as gen_interface
 import sharpy.utils.algebra as algebra
 import sharpy.linear.src.linuvlm as linuvlm
@@ -113,7 +113,7 @@ class StepLinearUVLM(BaseSolver):
     settings_default['cfl1'] = True
     settings_description['cfl1'] = 'If it is ``True``, it assumes that the discretisation complies with CFL=1'
     
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     scaling_settings_types = dict()
@@ -171,8 +171,8 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default, no_ctype=True)
-        su.to_custom_types(self.settings['ScalingDict'], self.scaling_settings_types,
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default, no_ctype=True)
+        settings_utils.to_custom_types(self.settings['ScalingDict'], self.scaling_settings_types,
                                  self.scaling_settings_default, no_ctype=True)
 
         # Initialise velocity generator
@@ -322,12 +322,12 @@ def run(self, **kwargs):
 
         """
 
-        aero_tstep = su.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
-        structure_tstep = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
-        convect_wake = su.set_value_or_default(kwargs, 'convect_wake', False)
-        dt= su.set_value_or_default(kwargs, 'dt', self.settings['dt'])                                                                                                    
-        t = su.set_value_or_default(kwargs, 't', self.data.ts*dt)
-        unsteady_contribution = su.set_value_or_default(kwargs, 'unsteady_contribution', False)
+        aero_tstep = settings_utils.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
+        structure_tstep = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        convect_wake = settings_utils.set_value_or_default(kwargs, 'convect_wake', False)
+        dt= settings_utils.set_value_or_default(kwargs, 'dt', self.settings['dt'])                                                                                                    
+        t = settings_utils.set_value_or_default(kwargs, 't', self.data.ts*dt)
+        unsteady_contribution = settings_utils.set_value_or_default(kwargs, 'unsteady_contribution', False)
 
         integr_order = self.settings['integr_order']
 
diff --git a/sharpy/solvers/stepuvlm.py b/sharpy/solvers/stepuvlm.py
index c17e3cf2b..d4cd9ad40 100644
--- a/sharpy/solvers/stepuvlm.py
+++ b/sharpy/solvers/stepuvlm.py
@@ -1,11 +1,9 @@
-import ctypes as ct
 import numpy as np
 import scipy.optimize
 import scipy.signal
 
-import sharpy.utils.algebra as algebra
 import sharpy.aero.utils.uvlmlib as uvlmlib
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.solver_interface import solver, BaseSolver
 import sharpy.utils.generator_interface as gen_interface
 import sharpy.utils.cout_utils as cout
@@ -133,7 +131,7 @@ class StepUvlm(BaseSolver):
     settings_default['quasi_steady'] = False
     settings_description['quasi_steady'] = 'Use quasi-steady approximation in UVLM'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
 
     def __init__(self):
@@ -150,7 +148,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            self.settings_options)
@@ -190,12 +188,12 @@ def run(self, **kwargs):
         """
 
         # Default values
-        aero_tstep = su.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
-        structure_tstep = su.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
-        convect_wake = su.set_value_or_default(kwargs, 'convect_wake', True)
-        dt= su.set_value_or_default(kwargs, 'dt', self.settings['dt'])
-        t = su.set_value_or_default(kwargs, 't', self.data.ts*dt)
-        unsteady_contribution = su.set_value_or_default(kwargs, 'unsteady_contribution', False)
+        aero_tstep = settings_utils.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
+        structure_tstep = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
+        convect_wake = settings_utils.set_value_or_default(kwargs, 'convect_wake', True)
+        dt= settings_utils.set_value_or_default(kwargs, 'dt', self.settings['dt'])
+        t = settings_utils.set_value_or_default(kwargs, 't', self.data.ts*dt)
+        unsteady_contribution = settings_utils.set_value_or_default(kwargs, 'unsteady_contribution', False)
 
         if not aero_tstep.zeta:
             return self.data
diff --git a/sharpy/solvers/timeintegrators.py b/sharpy/solvers/timeintegrators.py
index c60d08cc3..a11a96e5e 100644
--- a/sharpy/solvers/timeintegrators.py
+++ b/sharpy/solvers/timeintegrators.py
@@ -1,7 +1,7 @@
 import numpy as np
 import ctypes as ct
 
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.solver_interface import solver
 
 
@@ -81,7 +81,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            no_ctype=True)
@@ -180,7 +180,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                                  self.settings_types,
                                  self.settings_default,
                                  no_ctype=True)
diff --git a/sharpy/solvers/trim.py b/sharpy/solvers/trim.py
index af8dcaaeb..b8543f21e 100644
--- a/sharpy/solvers/trim.py
+++ b/sharpy/solvers/trim.py
@@ -1,13 +1,10 @@
-import ctypes as ct
 import numpy as np
 import scipy.optimize
-import warnings
 
-import sharpy.aero.utils.mapping as mapping
 import sharpy.utils.cout_utils as cout
 import sharpy.utils.solver_interface as solver_interface
 from sharpy.utils.solver_interface import solver, BaseSolver
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 import sharpy.utils.algebra as algebra
 
 
@@ -86,7 +83,7 @@ class Trim(BaseSolver):
     settings_default['refine_solution'] = False
     settings_description['refine_solution'] = 'If ``True`` and the optimiser routine allows for it, the optimiser will try to improve the solution with hybrid methods'
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -102,7 +99,7 @@ def __init__(self):
     def initialise(self, data, restart=False):
         self.data = data
         self.settings = data.settings[self.solver_id]
-        su.to_custom_types(self.settings, self.settings_types, self.settings_default)
+        settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
 
         self.solver = solver_interface.initialise_solver(self.settings['solver'])
         self.solver.initialise(self.data, self.settings['solver_settings'], restart=restart)
diff --git a/sharpy/solvers/updatepickle.py b/sharpy/solvers/updatepickle.py
index a3ca40f04..e23850cc5 100644
--- a/sharpy/solvers/updatepickle.py
+++ b/sharpy/solvers/updatepickle.py
@@ -1,5 +1,5 @@
 import numpy as np
-import sharpy.utils.settings as su
+import sharpy.utils.settings as settings_utils
 from sharpy.utils.solver_interface import solver, BaseSolver
 
 
@@ -14,7 +14,7 @@ class UpdatePickle(BaseSolver):
     settings_default = dict()
     settings_description = dict()
 
-    settings_table = su.SettingsTable()
+    settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
     def __init__(self):
@@ -27,7 +27,7 @@ def initialise(self, data, custom_settings=None, restart=False):
             self.settings = data.settings[self.solver_id]
         else:
             self.settings = custom_settings
-        su.to_custom_types(self.settings,
+        settings_utils.to_custom_types(self.settings,
                            self.settings_types,
                            self.settings_default,
                            no_ctype=True)
diff --git a/sharpy/utils/exceptions.py b/sharpy/utils/exceptions.py
index 8e889aef1..59a55c061 100644
--- a/sharpy/utils/exceptions.py
+++ b/sharpy/utils/exceptions.py
@@ -70,6 +70,17 @@ def __init__(self, setting, variable, options, value=None, message=''):
         message = 'The setting %s with entry %s is not one of the valid options: %s' % (setting, variable, options)
         super().__init__(variable, value, message=message)
         self.output_message(message, color_id = 4)
+        
+    
+class NotValidSettingType(DefaultValueBaseException):
+    """
+    Raised when a user gives a setting with an invalid type
+    """
+
+    def __init__(self, setting, variable,  data_types, value=None, message=''):
+        message = 'The setting %s with entry %s is not one of the valid types: %s' % (setting, variable, data_types)
+        super().__init__(variable, value, message=message)
+        self.output_message(message, color_id = 4)
 
 
 class SolverNotFound(Exception):
diff --git a/sharpy/utils/generate_cases.py b/sharpy/utils/generate_cases.py
index b2b5c74c1..e2b489fe2 100644
--- a/sharpy/utils/generate_cases.py
+++ b/sharpy/utils/generate_cases.py
@@ -252,7 +252,7 @@ def get_ielem_inode(connectivities, inode):
 
 def get_aoacl0_from_camber(x, y):
     """
-    This section provies the angle of attach of zero lift for a thin
+    This section provies the angle of attack of zero lift for a thin
     airfoil which camber line is defined by 'x' and 'y' coordinates
 
     Check Theory of wing sections. Abbott. pg 69
diff --git a/sharpy/utils/settings.py b/sharpy/utils/settings.py
index cdab4a4da..1e400e6b8 100644
--- a/sharpy/utils/settings.py
+++ b/sharpy/utils/settings.py
@@ -35,165 +35,16 @@ def cast(k, v, pytype, ctype, default):
 
 def to_custom_types(dictionary, types, default, options=dict(), no_ctype=True):
     for k, v in types.items():
-        if v == 'int':
-            if no_ctype:
-                data_type = int
-            else:
-                data_type = ct.c_int
-            try:
-                dictionary[k] = cast(k, dictionary[k], int, data_type, default[k])
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = cast(k, default[k], int, data_type, default[k])
-                notify_default_value(k, dictionary[k])
-
-        elif v == 'float':
-            if no_ctype:
-                data_type = float
-            else:
-                data_type = ct.c_double
-            try:
-                dictionary[k] = cast(k, dictionary[k], float, data_type, default[k])
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = cast(k, default[k], float, data_type, default[k])
-                notify_default_value(k, dictionary[k])
-
-        elif v == 'str':
-            try:
-                dictionary[k] = cast(k, dictionary[k], str, str, default[k])
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = cast(k, default[k], eval(v), eval(v), default[k])
-                notify_default_value(k, dictionary[k])
-
-        elif v == 'bool':
-            if no_ctype:
-                data_type = bool
-            else:
-                data_type = ct.c_bool
-            try:
-                dictionary[k] = cast(k, dictionary[k], str2bool, data_type, default[k])
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = cast(k, default[k], str2bool, data_type, default[k])
-                notify_default_value(k, dictionary[k])
-
-        elif v == 'list(str)':
-            try:
-                # if isinstance(dictionary[k], list):
-                #     continue
-                # dictionary[k] = dictionary[k].split(',')
-                # getting rid of leading and trailing spaces
-                dictionary[k] = list(map(lambda x: x.strip(), dictionary[k]))
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = default[k].copy()
-                notify_default_value(k, dictionary[k])
-
-        elif v == 'list(dict)':
-            try:
-                # if isinstance(dictionary[k], list):
-                #     continue
-                # dictionary[k] = dictionary[k].split(',')
-                # getting rid of leading and trailing spaces
-                for i in range(len(dictionary[k])):
-                    dictionary[k][i] = ast.literal_eval(dictionary[k][i])
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = default[k].copy()
-                notify_default_value(k, dictionary[k])
-
-        elif v == 'list(float)':
-            try:
-                dictionary[k]
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = default[k].copy()
-                notify_default_value(k, dictionary[k])
-
-            if isinstance(dictionary[k], np.ndarray):
-                continue
-            if isinstance(dictionary[k], list):
-                for i in range(len(dictionary[k])):
-                    dictionary[k][i] = float(dictionary[k][i])
-                dictionary[k] = np.array(dictionary[k])
-                continue
-            # dictionary[k] = dictionary[k].split(',')
-            # # getting rid of leading and trailing spaces
-            # dictionary[k] = list(map(lambda x: x.strip(), dictionary[k]))
-            if dictionary[k].find(',') < 0:
-                dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=' ', dtype=ct.c_double)
-            else:
-                dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=',', dtype=ct.c_double)
-
-        elif v == 'list(int)':
-            try:
-                dictionary[k]
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = default[k].copy()
-                notify_default_value(k, dictionary[k])
-
-            if isinstance(dictionary[k], np.ndarray):
-                continue
-            if isinstance(dictionary[k], list):
-                for i in range(len(dictionary[k])):
-                    dictionary[k][i] = int(dictionary[k][i])
-                dictionary[k] = np.array(dictionary[k])
-                continue
-            # dictionary[k] = dictionary[k].split(',')
-            # # getting rid of leading and trailing spaces
-            # dictionary[k] = list(map(lambda x: x.strip(), dictionary[k]))
-            if dictionary[k].find(',') < 0:
-                dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=' ').astype(ct.c_int)
-            else:
-                dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=',').astype(ct.c_int)
-
-        elif v == 'list(complex)':
-            try:
-                dictionary[k]
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = default[k].copy()
-                notify_default_value(k, dictionary[k])
-
-            if isinstance(dictionary[k], np.ndarray):
-                continue
-            if isinstance(dictionary[k], list):
-                for i in range(len(dictionary[k])):
-                    dictionary[k][i] = complex(dictionary[k][i])
-                dictionary[k] = np.array(dictionary[k])
-                continue
-            # dictionary[k] = dictionary[k].split(',')
-            # # getting rid of leading and trailing spaces
-            # dictionary[k] = list(map(lambda x: x.strip(), dictionary[k]))
-            if dictionary[k].find(',') < 0:
-                dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=' ').astype(complex)
-            else:
-                dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=',').astype(complex)
-
-        elif v == 'dict':
-            try:
-                if not isinstance(dictionary[k], dict):
-                    raise TypeError('Setting for {:s} is not a dictionary'.format(k))
-            except KeyError:
-                if default[k] is None:
-                    raise exceptions.NoDefaultValueException(k)
-                dictionary[k] = default[k].copy()
-                notify_default_value(k, dictionary[k])
+        if type(v) != list:
+            data_type = v
         else:
-            raise TypeError('Variable %s has an unknown type (%s) that cannot be casted' % (k, v))
-
+            if k in dictionary:
+                data_type = get_data_type_for_several_options(dictionary[k], v, k)
+            else:
+                # Choose first data type  in list for default value
+                data_type = v[0]
+        dictionary[k] = get_custom_type(dictionary, data_type, k, default, no_ctype)
+      
     check_settings_in_options(dictionary, types, options)
 
     unrecognised_settings = []
@@ -208,6 +59,168 @@ def to_custom_types(dictionary, types, default, options=dict(), no_ctype=True):
         raise Exception(unrecognised_settings)
 
 
+def get_data_type_for_several_options(dict_value, list_settings_types, setting_name):
+    """
+    Checks the data type of the setting input in case of several data type options. 
+    Only a scalar or list can be the case for these cases.  
+
+    Args:
+        dict_values: Dictionary value of processed settings
+        list_settings_types (list): Possible setting type options for this setting
+
+    Raises:
+        exception.NotValidSetting: if the setting is not allowed.
+    """
+    for data_type in list_settings_types:
+        if 'list' in data_type and (type(dict_value) == list or not np.isscalar(dict_value)):
+                return data_type
+        elif 'list' not in data_type and np.isscalar(dict_value):
+                return data_type
+    exceptions.NotValidSettingType(setting_name, dict_value, list_settings_types)
+
+def get_default_value(default_value, k, v, data_type = None, py_type = None):
+    if default_value is None:
+        raise exceptions.NoDefaultValueException(k)
+    if v in ['float', 'int', 'bool']:
+        converted_value = cast(k, default_value, py_type, data_type, default_value)
+    elif v == 'str':
+        converted_value = cast(k, default_value, eval(v), eval(v), default_value)
+    else:
+        converted_value = default_value.copy()
+    notify_default_value(k, converted_value)
+    return converted_value
+
+def get_custom_type(dictionary, v, k, default, no_ctype):
+    if v == 'int':
+        if no_ctype:
+            data_type = int
+        else:
+            data_type = ct.c_int
+        try:
+            dictionary[k] = cast(k, dictionary[k], int, data_type, default[k])
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v, data_type=data_type, py_type=int)
+
+    elif v == 'float':
+        if no_ctype:
+            data_type = float
+        else:
+            data_type = ct.c_double
+        try:
+            dictionary[k] = cast(k, dictionary[k], float, data_type, default[k])
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v, data_type=data_type, py_type=float)
+
+    elif v == 'str':
+        try:
+            dictionary[k] = cast(k, dictionary[k], str, str, default[k])
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v)
+
+    elif v == 'bool':
+        if no_ctype:
+            data_type = bool
+        else:
+            data_type = ct.c_bool
+        try:
+            dictionary[k] = cast(k, dictionary[k], str2bool, data_type, default[k])
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v, data_type=data_type, py_type=str2bool)
+
+    elif v == 'list(str)':
+        try:
+            # if isinstance(dictionary[k], list):
+            #     continue
+            # dictionary[k] = dictionary[k].split(',')
+            # getting rid of leading and trailing spaces
+            dictionary[k] = list(map(lambda x: x.strip(), dictionary[k]))
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v)
+
+    elif v == 'list(dict)':
+        try:
+            # if isinstance(dictionary[k], list):
+            #     continue
+            # dictionary[k] = dictionary[k].split(',')
+            # getting rid of leading and trailing spaces
+            for i in range(len(dictionary[k])):
+                dictionary[k][i] = ast.literal_eval(dictionary[k][i])
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v)
+
+    elif v == 'list(float)':
+        try:
+            dictionary[k]
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v)
+
+        if isinstance(dictionary[k], np.ndarray):
+            return dictionary[k]
+        if isinstance(dictionary[k], list):
+            for i in range(len(dictionary[k])):
+                dictionary[k][i] = float(dictionary[k][i])
+            dictionary[k] = np.array(dictionary[k])
+            return dictionary[k]
+        # dictionary[k] = dictionary[k].split(',')
+        # # getting rid of leading and trailing spaces
+        # dictionary[k] = list(map(lambda x: x.strip(), dictionary[k]))
+        if dictionary[k].find(',') < 0:
+            dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=' ', dtype=ct.c_double)
+        else:
+            dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=',', dtype=ct.c_double)
+
+    elif v == 'list(int)':
+        try:
+            dictionary[k]
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v)
+
+        if isinstance(dictionary[k], np.ndarray):
+            return dictionary[k]
+        if isinstance(dictionary[k], list):
+            for i in range(len(dictionary[k])):
+                dictionary[k][i] = int(dictionary[k][i])
+            dictionary[k] = np.array(dictionary[k])
+            return dictionary[k]
+        # dictionary[k] = dictionary[k].split(',')
+        # # getting rid of leading and trailing spaces
+        # dictionary[k] = list(map(lambda x: x.strip(), dictionary[k]))
+        if dictionary[k].find(',') < 0:
+            dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=' ').astype(ct.c_int)
+        else:
+            dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=',').astype(ct.c_int)
+
+    elif v == 'list(complex)':
+        try:
+            dictionary[k]
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v)
+
+        if isinstance(dictionary[k], np.ndarray):
+            return dictionary[k]
+        if isinstance(dictionary[k], list):
+            for i in range(len(dictionary[k])):
+                dictionary[k][i] = complex(dictionary[k][i])
+            dictionary[k] = np.array(dictionary[k])
+            return dictionary[k]
+        # dictionary[k] = dictionary[k].split(',')
+        # # getting rid of leading and trailing spaces
+        # dictionary[k] = list(map(lambda x: x.strip(), dictionary[k]))
+        if dictionary[k].find(',') < 0:
+            dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=' ').astype(complex)
+        else:
+            dictionary[k] = np.fromstring(dictionary[k].strip('[]'), sep=',').astype(complex)
+
+    elif v == 'dict':
+        try:
+            if not isinstance(dictionary[k], dict):
+                raise TypeError('Setting for {:s} is not a dictionary'.format(k))
+        except KeyError:
+            dictionary[k] = get_default_value(default[k], k, v)
+    else:
+        raise TypeError('Variable %s has an unknown type (%s) that cannot be casted' % (k, v))
+    return dictionary[k]
+
 def check_settings_in_options(settings, settings_types, settings_options):
     """
     Checks that settings given a type ``str`` or ``int`` and allowable options are indeed valid.


Test patch:
diff --git a/tests/coupled/static/pazy/generate_pazy.py b/tests/coupled/static/pazy/generate_pazy.py
index f756915d7..9dc1b475d 100644
--- a/tests/coupled/static/pazy/generate_pazy.py
+++ b/tests/coupled/static/pazy/generate_pazy.py
@@ -13,6 +13,7 @@ def generate_pazy(u_inf, case_name, output_folder='/output/', cases_folder='', *
     M = kwargs.get('M', 4)
     N = kwargs.get('N', 32)
     M_star_fact = kwargs.get('Ms', 10)
+    flag_multiple_mstar_input = kwargs.get('test_multiple_inputs', False)
 
     # SHARPy nonlinear reference solution
     ws = wings.PazyControlSurface(M=M,
@@ -66,64 +67,67 @@ def generate_pazy(u_inf, case_name, output_folder='/output/', cases_folder='', *
         'wake_shape_generator_input': {'u_inf': ws.u_inf,
                                        'u_inf_direction': ws.u_inf_direction,
                                        'dt': ws.dt}}
-
-    ws.config['StaticUvlm'] = {
-        'rho': ws.rho,
-        'velocity_field_generator': 'SteadyVelocityField',
-        'velocity_field_input': {
-            'u_inf': ws.u_inf,
-            'u_inf_direction': ws.u_inf_direction},
-        'rollup_dt': ws.dt,
-        'print_info': 'on',
-        'horseshoe': 'on',
-        'num_cores': 4,
-        'n_rollup': 0,
-        'rollup_aic_refresh': 0,
-        'rollup_tolerance': 1e-4}
-
-    settings = dict()
-    settings['NonLinearStatic'] = {'print_info': 'off',
-                                   'max_iterations': 200,
-                                   'num_load_steps': 5,
-                                   'delta_curved': 1e-6,
-                                   'min_delta': 1e-8,
-                                   'gravity_on': gravity_on,
-                                   'gravity': 9.81}
-
-    ws.config['StaticCoupled'] = {
-        'print_info': 'on',
-        'max_iter': 200,
-        'n_load_steps': 4,
-        'tolerance': 1e-5,
-        'relaxation_factor': 0.1,
-        'aero_solver': 'StaticUvlm',
-        'aero_solver_settings': {
+    if flag_multiple_mstar_input:
+        ws.config['AerogridLoader']['mstar'] = [ws.config['AerogridLoader']['mstar'], ws.config['AerogridLoader']['mstar']]
+        ws.config['SHARPy']['flow'] = ['BeamLoader','AerogridLoader']
+    else:
+        ws.config['StaticUvlm'] = {
             'rho': ws.rho,
-            'print_info': 'off',
-            'horseshoe': 'on',
-            'num_cores': 4,
-            'n_rollup': 0,
-            'rollup_dt': ws.dt,
-            'rollup_aic_refresh': 1,
-            'rollup_tolerance': 1e-4,
             'velocity_field_generator': 'SteadyVelocityField',
             'velocity_field_input': {
                 'u_inf': ws.u_inf,
                 'u_inf_direction': ws.u_inf_direction},
-            'vortex_radius': 1e-9},
-        'structural_solver': 'NonLinearStatic',
-        'structural_solver_settings': settings['NonLinearStatic']}
-
-    ws.config['AerogridPlot'] = {'include_rbm': 'off',
-                                 'include_applied_forces': 'on',
-                                 'minus_m_star': 0}
-
-    ws.config['BeamPlot'] = {'include_rbm': 'off',
-                             'include_applied_forces': 'on'}
-
-    ws.config['WriteVariablesTime'] = {'structure_variables': ['pos'],
-                                       'structure_nodes': list(range(0, ws.num_node_surf)),
-                                       'cleanup_old_solution': 'on'}
+            'rollup_dt': ws.dt,
+            'print_info': 'on',
+            'horseshoe': 'on',
+            'num_cores': 4,
+            'n_rollup': 0,
+            'rollup_aic_refresh': 0,
+            'rollup_tolerance': 1e-4}
+
+        settings = dict()
+        settings['NonLinearStatic'] = {'print_info': 'off',
+                                    'max_iterations': 200,
+                                    'num_load_steps': 5,
+                                    'delta_curved': 1e-6,
+                                    'min_delta': 1e-8,
+                                    'gravity_on': gravity_on,
+                                    'gravity': 9.81}
+
+        ws.config['StaticCoupled'] = {
+            'print_info': 'on',
+            'max_iter': 200,
+            'n_load_steps': 4,
+            'tolerance': 1e-5,
+            'relaxation_factor': 0.1,
+            'aero_solver': 'StaticUvlm',
+            'aero_solver_settings': {
+                'rho': ws.rho,
+                'print_info': 'off',
+                'horseshoe': 'on',
+                'num_cores': 4,
+                'n_rollup': 0,
+                'rollup_dt': ws.dt,
+                'rollup_aic_refresh': 1,
+                'rollup_tolerance': 1e-4,
+                'velocity_field_generator': 'SteadyVelocityField',
+                'velocity_field_input': {
+                    'u_inf': ws.u_inf,
+                    'u_inf_direction': ws.u_inf_direction},
+                'vortex_radius': 1e-9},
+            'structural_solver': 'NonLinearStatic',
+            'structural_solver_settings': settings['NonLinearStatic']}
+
+        ws.config['AerogridPlot'] = {'include_rbm': 'off',
+                                    'include_applied_forces': 'on',
+                                    'minus_m_star': 0}
+
+        ws.config['BeamPlot'] = {'include_rbm': 'off',
+                                'include_applied_forces': 'on'}
+
+        ws.config['WriteVariablesTime'] = {'structure_variables': ['pos'],
+                                        'structure_nodes': list(range(0, ws.num_node_surf)),
+                                        'cleanup_old_solution': 'on'}
 
     ws.config.write()
 
diff --git a/tests/utils/test_settings.py b/tests/utils/test_settings.py
index aad222495..2a2bac4a7 100644
--- a/tests/utils/test_settings.py
+++ b/tests/utils/test_settings.py
@@ -4,6 +4,8 @@
 from copy import deepcopy
 import numpy as np
 import unittest
+import tests.coupled.static.pazy.generate_pazy as gp
+import os
 
 
 class TestSettings(unittest.TestCase):
@@ -11,6 +13,8 @@ class TestSettings(unittest.TestCase):
     Tests the settings utilities module
     """
 
+    route_test_dir = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))
+    
     def setUp(self):
         cout.start_writer()
 
@@ -146,3 +150,29 @@ def test_settings_to_custom_types(self):
             temp_in_dict = deepcopy(original_dict)
             del temp_in_dict['float_list_var']
             result = settings.to_custom_types(temp_in_dict, types_dict, temp_default_dict)
+
+    def test_multiple_setting_option(self):
+        u_inf = 50
+        alpha = 7
+        case_name = 'pazy_uinf{:04g}_alpha{:04g}'.format(u_inf * 10, alpha * 10)
+
+        M = 16
+        N = 64
+        Msf = 1
+
+        cases_folder = self.route_test_dir + '/pazy/cases/'
+        output_folder = self.route_test_dir + '/pazy/cases/'
+        # run case
+        gp.generate_pazy(u_inf, case_name, output_folder, cases_folder,
+                         alpha=alpha,
+                         M=M,
+                         N=N,
+                         Msf=Msf,
+                         test_multiple_inputs=True)
+
+    def tearDown(self):
+        cases_folder = self.route_test_dir + '/pazy/cases/'
+
+        if os.path.isdir(cases_folder):
+            import shutil
+            shutil.rmtree(cases_folder)