Repository: ImperialCollegeLondon/sharpy (language: python)

Issue:
Include gravity direction as input for structural solvers
Right now gravity always acts on the z direction in SHARPy.

However, xbeam allows gravity to be applied in any direction. Thus, it seems that we have ready the most laborious part of this feature and we are just missing the interface. 

New feature: apply external forces to the system at runtime
This new feature allow the application of external forces to the system computed at runtime based on the state of the system at each time step. 

In particular, forces related to floating dynamics of wind turbines (mooring lines, hydrostatics, hydrodynamics and wave forces) are considered.

This PR Fixes #112 including the possibility to change the direction of gravity.

It also includes other minor changes:
- Include lumped masses defined by a 6x6 mass matrix through the structural input variable ``lumped_mass_mat``
- Possibility to skip the writing of the aerodynamic wake in ``SaveData`` through the ``save_wake`` setting.
- Generate spar floating wind turbines with the ``template_wt`` library.

It updates some packages in the conda environment. 

Requirements before merging
- [x] Update UVLM library
- [x] Clean up
- [x] Merge PR #125 (multibody) before this one



Discussion / hints:
Now the gravity direction is a property of the simulation and not a beam property.

It is defined through the "gravity_dir" setting in structural solvers.

Everything should be backwards compatible except from solvers calling "cbeam3_asbly_static". I think the call to this function was not right. To my knowledge, this only affect static multibody computations. 

I do not run many computations with gravity so any test or feedback is welcome. 


Fix patch:
diff --git a/cases/templates/template_wt.py b/cases/templates/template_wt.py
index e18db01d6..6fe674fee 100644
--- a/cases/templates/template_wt.py
+++ b/cases/templates/template_wt.py
@@ -88,6 +88,212 @@ def create_blade_coordinates(num_node, node_r, node_y, node_z):
     return coordinates
 
 
+######################################################################
+# FROM excel type04
+######################################################################
+def spar_from_excel_type04(op_params,
+                            geom_params,
+                            excel_description,
+                            options):
+    """
+    spar_from_excel_type04
+
+    Function needed to generate a spar floating wind turbine rotor from an excel database type04.
+    Rotor + tower + floating spar
+
+    Args:
+        op_param (dict): Dictionary with operating parameters
+        geom_param (dict): Dictionray with geometical parameters
+        excel_description (dict): Dictionary describing the sheets of the excel file
+        option (dict): Dictionary with the different options for the wind turbine generation
+
+    Returns:
+        floating (sharpy.utils.generate_cases.AeroelasticInfromation): Aeroelastic information of the spar floating wind turbine
+    """
+
+    # Generate the tower + rotor
+    wt, LC, MB, hub_node = generate_from_excel_type03(op_params,
+                                    geom_params,
+                                    excel_description,
+                                    options)
+
+    # Remove base clam
+    wt.StructuralInformation.boundary_conditions[0] = 0
+
+    excel_file_name = excel_description['excel_file_name']
+    excel_sheet_parameters = excel_description['excel_sheet_parameters']
+    excel_sheet_structural_spar = excel_description['excel_sheet_structural_spar']
+    tol_remove_points = geom_params['tol_remove_points']
+
+    # Generate the spar
+    SparHeight = gc.read_column_sheet_type01(excel_file_name,
+                                             excel_sheet_parameters,
+                                             'SparHeight')
+    SparBallastMass = gc.read_column_sheet_type01(excel_file_name,
+                                                  excel_sheet_parameters,
+                                                  'SparBallastMass')
+    SparBallastDepth = gc.read_column_sheet_type01(excel_file_name,
+                                                   excel_sheet_parameters,
+                                                   'SparBallastDepth')
+
+    # Read from excel file
+    SparHtFract = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparHtFract')
+    SparMassDen = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparMassDen')
+    SparFAStif = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparFAStif')
+    SparSSStif = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparSSStif')
+    SparGJStif = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparGJStif')
+    SparEAStif = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparEAStif')
+    SparFAIner = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparFAIner')
+    SparSSIner = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparSSIner')
+    SparFAcgOf = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparFAcgOf')
+    SparSScgOf = gc.read_column_sheet_type01(excel_file_name, excel_sheet_structural_spar, 'SparSScgOf')
+
+    ElevationSpar = SparHtFract*SparHeight
+    spar = gc.AeroelasticInformation()
+    spar.StructuralInformation.num_elem = len(ElevationSpar) - 2
+    spar.StructuralInformation.num_node_elem = 3
+    spar.StructuralInformation.compute_basic_num_node()
+
+    # Interpolate excel variables into the correct locations
+    node_r, elem_r = create_node_radial_pos_from_elem_centres(ElevationSpar,
+                                    spar.StructuralInformation.num_node,
+                                    spar.StructuralInformation.num_elem,
+                                    spar.StructuralInformation.num_node_elem)
+
+    # Stiffness
+    elem_EA = np.interp(elem_r, ElevationSpar, SparEAStif)
+    elem_EIz = np.interp(elem_r, ElevationSpar, SparSSStif)
+    elem_EIy = np.interp(elem_r, ElevationSpar, SparFAStif)
+    elem_GJ = np.interp(elem_r, ElevationSpar, SparGJStif)
+    # Stiffness: estimate unknown properties
+    cout.cout_wrap.print_file = False
+    cout.cout_wrap('WARNING: The poisson cofficient is assumed equal to 0.3', 3)
+    cout.cout_wrap('WARNING: Cross-section area is used as shear area', 3)
+    poisson_coef = 0.3
+    elem_GAy = elem_EA/2.0/(1.0+poisson_coef)
+    elem_GAz = elem_EA/2.0/(1.0+poisson_coef)
+
+    # Inertia
+    elem_mass_per_unit_length = np.interp(elem_r, ElevationSpar, SparMassDen)
+    elem_mass_iner_y = np.interp(elem_r, ElevationSpar, SparFAIner)
+    elem_mass_iner_z = np.interp(elem_r, ElevationSpar, SparSSIner)
+    # TODO: check yz axis and Flap-edge
+    elem_pos_cg_B = np.zeros((spar.StructuralInformation.num_elem, 3),)
+    elem_pos_cg_B[:, 1] = np.interp(elem_r, ElevationSpar, SparSScgOf)
+    elem_pos_cg_B[:, 2] = np.interp(elem_r, ElevationSpar, SparFAcgOf)
+
+    # Stiffness: estimate unknown properties
+    cout.cout_wrap('WARNING: Using perpendicular axis theorem to compute the inertia around xB', 3)
+    elem_mass_iner_x = elem_mass_iner_y + elem_mass_iner_z
+
+    # Create the tower
+    spar.StructuralInformation.create_mass_db_from_vector(elem_mass_per_unit_length, elem_mass_iner_x, elem_mass_iner_y, elem_mass_iner_z, elem_pos_cg_B)
+    spar.StructuralInformation.create_stiff_db_from_vector(elem_EA, elem_GAy, elem_GAz, elem_GJ, elem_EIy, elem_EIz)
+
+    coordinates = np.zeros((spar.StructuralInformation.num_node, 3),)
+    coordinates[:, 0] = node_r
+
+    spar.StructuralInformation.generate_1to1_from_vectors(
+        num_node_elem=spar.StructuralInformation.num_node_elem,
+        num_node=spar.StructuralInformation.num_node,
+        num_elem=spar.StructuralInformation.num_elem,
+        coordinates=coordinates,
+        stiffness_db=spar.StructuralInformation.stiffness_db,
+        mass_db=spar.StructuralInformation.mass_db,
+        frame_of_reference_delta='y_AFoR',
+        vec_node_structural_twist=np.zeros((spar.StructuralInformation.num_node,),),
+        num_lumped_mass=1)
+
+    spar.StructuralInformation.boundary_conditions[0] = -1
+    spar.StructuralInformation.boundary_conditions[-1] = 1
+
+    # Include ballast mass
+    dist = np.abs(coordinates[:, 0] + SparBallastDepth)
+    min_dist = np.amin(dist)
+    loc_min = np.where(dist == min_dist)[0]
+
+    cout.cout_wrap("Ballast at node %d at position %f" % (loc_min, coordinates[loc_min, 0]), 2)
+
+    spar.StructuralInformation.lumped_mass_nodes[0] = loc_min
+    spar.StructuralInformation.lumped_mass[0] = SparBallastMass
+    spar.StructuralInformation.lumped_mass_inertia[0] = np.zeros((3,3))
+    spar.StructuralInformation.lumped_mass_position[0] = np.zeros((3,))
+
+    spar.AerodynamicInformation.set_to_zero(spar.StructuralInformation.num_node_elem,
+                                        spar.StructuralInformation.num_node,
+                                        spar.StructuralInformation.num_elem)
+
+    # Generate tower base
+    TowerBaseHeight = gc.read_column_sheet_type01(excel_file_name,
+                                                  excel_sheet_parameters,
+                                                  'TowerBaseHeight')
+    NodesBase = gc.read_column_sheet_type01(excel_file_name,
+                                                  excel_sheet_parameters,
+                                                  'NodesBase')
+
+
+    base = gc.AeroelasticInformation()
+    base.StructuralInformation.num_node = NodesBase
+    base.StructuralInformation.num_node_elem = 3
+    base.StructuralInformation.compute_basic_num_elem()
+
+    node_coord = np.zeros((base.StructuralInformation.num_node, 3))
+    node_coord[:, 0] = np.linspace(0., TowerBaseHeight, base.StructuralInformation.num_node)
+
+    base_stiffness_db = np.zeros((base.StructuralInformation.num_elem, 6, 6))
+    base_mass_db = np.zeros((base.StructuralInformation.num_elem, 6, 6))
+    vec_node_structural_twist = np.zeros((base.StructuralInformation.num_node,))
+
+    for ielem in range(base.StructuralInformation.num_elem):
+        base_stiffness_db[ielem, :, :] = spar.StructuralInformation.stiffness_db[-1, :, :]
+        base_mass_db[ielem, :, :] = spar.StructuralInformation.mass_db[-1, :, :]
+    # for ielem in range(base.StructuralInformation.num_elem):
+    #     inode_cent = 2*ielem + 1
+    #     rel_dist_to_spar = node_coord[inode_cent, 0]/TowerBaseHeight
+    #     rel_dist_to_base = (node_coord[-1, 0] - node_coord[inode_cent, 0])/TowerBaseHeight
+    #     base_stiffness_db[ielem, :, :] = (spar.StructuralInformation.stiffness_db[0, :, :]*rel_dist_to_base +
+    #                                   wt.StructuralInformation.stiffness_db[0, :, :]*rel_dist_to_spar)
+    #     base_mass_db[ielem, :, :] = (spar.StructuralInformation.mass_db[0, :, :]*rel_dist_to_base +
+    #                                   wt.StructuralInformation.mass_db[0, :, :]*rel_dist_to_spar)
+
+    base.StructuralInformation.generate_1to1_from_vectors(
+                            base.StructuralInformation.num_node_elem,
+                            base.StructuralInformation.num_node,
+                            base.StructuralInformation.num_elem,
+                            node_coord,
+                            base_stiffness_db,
+                            base_mass_db,
+                            'y_AFoR',
+                            vec_node_structural_twist,
+                            num_lumped_mass=0)
+
+    base.StructuralInformation.boundary_conditions[0] = 1
+    base.StructuralInformation.body_number *= 0
+
+    base.AerodynamicInformation.set_to_zero(base.StructuralInformation.num_node_elem,
+                                        base.StructuralInformation.num_node,
+                                        base.StructuralInformation.num_elem)
+
+    # Assembly
+    spar.assembly(base)
+    spar.remove_duplicated_points(1e-6)
+    spar.StructuralInformation.body_number *= 0
+    nodes_spar = spar.StructuralInformation.num_node + 0
+    skip = [hub_node + nodes_spar]
+    wt.StructuralInformation.coordinates[:, 0] += TowerBaseHeight
+    spar.assembly(wt)
+    spar.remove_duplicated_points(1e-6, skip=skip)
+    for ielem in range(spar.StructuralInformation.num_elem):
+        if not spar.StructuralInformation.body_number[ielem] == 0:
+            spar.StructuralInformation.body_number[ielem] -= 1
+
+    # Update Lagrange Constraints and Multibody information
+    LC[0].node_in_body += nodes_spar - 1
+    MB[0].FoR_movement = 'free'
+    MB[1].FoR_position[0] += TowerBaseHeight
+
+    return spar, LC, MB
+
 ######################################################################
 # FROM excel type03
 ######################################################################
@@ -96,9 +302,9 @@ def rotor_from_excel_type03(in_op_params,
                             in_excel_description,
                             in_options):
     """
-    generate_from_excel_type03_db
+    rotor_from_excel_type03
 
-    Function needed to generate a wind turbine from an excel database type03
+    Function needed to generate a wind turbine rotor from an excel database type03
 
     Args:
         op_param (dict): Dictionary with operating parameters
@@ -495,70 +701,37 @@ def rotor_from_excel_type03(in_op_params,
     return rotor
 
 
-def generate_from_excel_type02(chord_panels,
-                                rotation_velocity,
-                                pitch_deg,
-                                excel_file_name='database_excel_type02.xlsx',
-                                excel_sheet_parameters='parameters',
-                                excel_sheet_structural_blade='structural_blade',
-                                excel_sheet_discretization_blade='discretization_blade',
-                                excel_sheet_aero_blade='aero_blade',
-                                excel_sheet_airfoil_info='airfoil_info',
-                                excel_sheet_airfoil_coord='airfoil_coord',
-                                excel_sheet_structural_tower='structural_tower',
-                                m_distribution='uniform',
-                                h5_cross_sec_prop=None,
-                                n_points_camber=100,
-                                tol_remove_points=1e-3,
-                                user_defined_m_distribution_type=None,
-                                wsp=0.,
-                                dt=0.):
+def generate_from_excel_type03(op_params,
+                                   geom_params,
+                                   excel_description,
+                                   options):
+
     """
-    generate_from_excel_type02
+    generate_from_excel_type03
 
-    Function needed to generate a wind turbine from an excel database according to OpenFAST inputs
+    Function needed to generate a wind turbine (tower + rotor) from an excel database type03.
+    See ``rotor_from_excel_type03'' for more information.
 
     Args:
-        chord_panels (int): Number of panels on the blade surface in the chord direction
-        rotation_velocity (float): Rotation velocity of the rotor
-        pitch_deg (float): pitch angle in degrees
-        excel_file_name (str):
-        excel_sheet_structural_blade (str):
-        excel_sheet_aero_blade (str):
-        excel_sheet_airfoil_coord (str):
-        excel_sheet_parameters (str):
-        excel_sheet_structural_tower (str):
-        m_distribution (str):
-        n_points_camber (int): number of points to define the camber of the airfoil,
-        tol_remove_points (float): maximum distance to remove adjacent points
-        user_defined_m_distribution_type (string): type of distribution of the chordwise panels when 'm_distribution' == 'user_defined'
-        camber_effects_on_twist (bool): When true plain airfoils are used and the blade is twisted and preloaded based on thin airfoil theory
-        wsp (float): wind speed (It may be needed for discretisation purposes)
-        dt (float): time step (It may be needed for discretisation purposes)
+        op_param (dict): Dictionary with operating parameters
+        geom_param (dict): Dictionray with geometical parameters
+        excel_description (dict): Dictionary describing the sheets of the excel file
+        option (dict): Dictionary with the different options for the wind turbine generation
 
     Returns:
-        wt (sharpy.utils.generate_cases.AeroelasticInfromation): Aeroelastic infrmation of the wind turbine
-        LC (list): list of all the Lagrange constraints needed in the cases (sharpy.utils.generate_cases.LagrangeConstraint)
-        MB (list): list of the multibody information of each body (sharpy.utils.generate_cases.BodyInfrmation)
+        wt (sharpy.utils.generate_cases.AeroelasticInfromation): Aeroelastic information of the wind turbine (tower + rotor)
     """
+    rotor = rotor_from_excel_type03(op_params,
+                                    geom_params,
+                                    excel_description,
+                                    options)
 
-    rotor = rotor_from_excel_type02(chord_panels,
-                                  rotation_velocity,
-                                  pitch_deg,
-                                  excel_file_name=excel_file_name,
-                                  excel_sheet_parameters=excel_sheet_parameters,
-                                  excel_sheet_structural_blade=excel_sheet_structural_blade,
-                                  excel_sheet_discretization_blade=excel_sheet_discretization_blade,
-                                  excel_sheet_aero_blade=excel_sheet_aero_blade,
-                                  excel_sheet_airfoil_info=excel_sheet_airfoil_info,
-                                  excel_sheet_airfoil_coord=excel_sheet_airfoil_coord,
-                                  m_distribution=m_distribution,
-                                  h5_cross_sec_prop=h5_cross_sec_prop,
-                                  n_points_camber=n_points_camber,
-                                  tol_remove_points=tol_remove_points,
-                                  user_defined_m_distribution_type=user_defined_m_distribution_type,
-                                  wsp=0.,
-                                  dt=0.)
+
+    excel_file_name = excel_description['excel_file_name']
+    excel_sheet_parameters = excel_description['excel_sheet_parameters']
+    excel_sheet_structural_tower = excel_description['excel_sheet_structural_tower']
+    tol_remove_points = geom_params['tol_remove_points']
+    rotation_velocity = op_params['rotation_velocity']
 
     ######################################################################
     ## TOWER
@@ -641,13 +814,19 @@ def generate_from_excel_type02(chord_panels,
 
     # Read overhang and nacelle properties from excel file
     overhang_len = gc.read_column_sheet_type01(excel_file_name, excel_sheet_parameters, 'overhang')
-    # HubMass = gc.read_column_sheet_type01(excel_file_name, excel_sheet_nacelle, 'HubMass')
     NacelleMass = gc.read_column_sheet_type01(excel_file_name, excel_sheet_parameters, 'NacMass')
+    NacelleMass_x = gc.read_column_sheet_type01(excel_file_name, excel_sheet_parameters, 'NacMass_x')
+    NacelleMass_z = gc.read_column_sheet_type01(excel_file_name, excel_sheet_parameters, 'NacMass_z')
     # NacelleYawIner = gc.read_column_sheet_type01(excel_file_name, excel_sheet_nacelle, 'NacelleYawIner')
 
     # Include nacelle mass
     tower.StructuralInformation.lumped_mass_nodes = np.array([tower.StructuralInformation.num_node - 1], dtype=int)
     tower.StructuralInformation.lumped_mass = np.array([NacelleMass], dtype=float)
+    if not NacelleMass_x is None and not NacelleMass_z is None:
+        tower.StructuralInformation.lumped_mass_position = np.array([np.array([NacelleMass_z, 0, NacelleMass_x])], dtype=float)
+    else:
+        cout.cout_wrap('WARNING: Nacelle mass placed at tower top', 3)
+
 
     tower.AerodynamicInformation.set_to_zero(tower.StructuralInformation.num_node_elem,
                                             tower.StructuralInformation.num_node,
@@ -657,6 +836,7 @@ def generate_from_excel_type02(chord_panels,
     # numberOfBlades = gc.read_column_sheet_type01(excel_file_name, excel_sheet_parameters, 'NumBl')
     tilt = gc.read_column_sheet_type01(excel_file_name, excel_sheet_parameters, 'ShftTilt')*deg2rad
     # cone = gc.read_column_sheet_type01(excel_file_name, excel_sheet_parameters, 'Cone')*deg2rad
+    HubMass = gc.read_column_sheet_type01(excel_file_name, excel_sheet_parameters, 'HubMass')
 
     overhang = gc.AeroelasticInformation()
     overhang.StructuralInformation.num_node = 3
@@ -664,8 +844,8 @@ def generate_from_excel_type02(chord_panels,
     overhang.StructuralInformation.compute_basic_num_elem()
     node_pos = np.zeros((overhang.StructuralInformation.num_node, 3), )
     node_pos[:, 0] += tower.StructuralInformation.coordinates[-1, 0]
-    node_pos[:, 0] += np.linspace(0., overhang_len*np.sin(tilt*deg2rad), overhang.StructuralInformation.num_node)
-    node_pos[:, 2] = np.linspace(0., -overhang_len*np.cos(tilt*deg2rad), overhang.StructuralInformation.num_node)
+    node_pos[:, 0] += np.linspace(0., -overhang_len*np.sin(tilt*deg2rad), overhang.StructuralInformation.num_node)
+    node_pos[:, 2] = np.linspace(0., overhang_len*np.cos(tilt*deg2rad), overhang.StructuralInformation.num_node)
     # TODO: change the following by real values
     # Same properties as the last element of the tower
     cout.cout_wrap("WARNING: Using the structural properties of the last tower section for the overhang", 3)
@@ -675,6 +855,12 @@ def generate_from_excel_type02(chord_panels,
     oh_GA = tower.StructuralInformation.stiffness_db[-1, 1, 1]
     oh_GJ = tower.StructuralInformation.stiffness_db[-1, 3, 3]
     oh_EI = tower.StructuralInformation.stiffness_db[-1, 4, 4]
+    if not HubMass is None:
+        num_lumped_mass_overhang = 1
+    else:
+        cout.cout_wrap('WARNING: HubMass not found', 3)
+        num_lumped_mass_overhang = 0
+
     overhang.StructuralInformation.generate_uniform_sym_beam(node_pos,
                                                             oh_mass_per_unit_length,
                                                             oh_mass_iner,
@@ -684,17 +870,24 @@ def generate_from_excel_type02(chord_panels,
                                                             oh_EI,
                                                             num_node_elem=3,
                                                             y_BFoR='y_AFoR',
-                                                            num_lumped_mass=0)
+                                                            num_lumped_mass=num_lumped_mass_overhang)
 
     overhang.StructuralInformation.boundary_conditions = np.zeros((overhang.StructuralInformation.num_node), dtype=int)
     overhang.StructuralInformation.boundary_conditions[-1] = -1
 
+    if not HubMass is None:
+        # Include hub mass
+        overhang.StructuralInformation.lumped_mass_nodes = np.array([overhang.StructuralInformation.num_node - 1], dtype=int)
+        overhang.StructuralInformation.lumped_mass = np.array([HubMass], dtype=float)
+
     overhang.AerodynamicInformation.set_to_zero(overhang.StructuralInformation.num_node_elem,
                                                 overhang.StructuralInformation.num_node,
                                                 overhang.StructuralInformation.num_elem)
 
     tower.assembly(overhang)
     tower.remove_duplicated_points(tol_remove_points)
+    hub_node = tower.StructuralInformation.num_node
+    tower.StructuralInformation.body_number *= 0
 
     ######################################################################
     ##  WIND TURBINE
@@ -747,7 +940,7 @@ def generate_from_excel_type02(chord_panels,
     ######################################################################
     ## RETURN
     ######################################################################
-    return wt, LC, MB
+    return wt, LC, MB, hub_node
 
 
 ######################################################################
@@ -810,3 +1003,63 @@ def rotor_from_excel_type02(chord_panels,
                                    options)
 
     return rotor
+
+
+def generate_from_excel_type02(chord_panels,
+                                rotation_velocity,
+                                pitch_deg,
+                                excel_file_name='database_excel_type02.xlsx',
+                                excel_sheet_parameters='parameters',
+                                excel_sheet_structural_blade='structural_blade',
+                                excel_sheet_discretization_blade='discretization_blade',
+                                excel_sheet_aero_blade='aero_blade',
+                                excel_sheet_airfoil_info='airfoil_info',
+                                excel_sheet_airfoil_coord='airfoil_coord',
+                                excel_sheet_structural_tower='structural_tower',
+                                m_distribution='uniform',
+                                h5_cross_sec_prop=None,
+                                n_points_camber=100,
+                                tol_remove_points=1e-3,
+                                user_defined_m_distribution_type=None,
+                                camber_effect_on_twist=False,
+                                wsp=0.,
+                                dt=0.):
+
+    # Warning for back compatibility
+    cout.cout_wrap('generate_from_excel_type02 is obsolete! rotor_from_excel_type03 instead!', 3)
+
+    # Assign values to dictionaries
+    op_params = {}
+    op_params['rotation_velocity'] = rotation_velocity
+    op_params['pitch_deg'] = pitch_deg
+    op_params['wsp'] = wsp
+    op_params['dt'] = dt
+
+    geom_params = {}
+    geom_params['chord_panels'] = chord_panels
+    geom_params['tol_remove_points'] = tol_remove_points
+    geom_params['n_points_camber'] = n_points_camber
+    geom_params['h5_cross_sec_prop'] = h5_cross_sec_prop
+    geom_params['m_distribution'] = m_distribution
+
+    options = {}
+    options['camber_effect_on_twist'] = camber_effect_on_twist
+    options['user_defined_m_distribution_type'] = user_defined_m_distribution_type
+    options['include_polars'] = False
+
+    excel_description = {}
+    excel_description['excel_file_name'] = excel_file_name
+    excel_description['excel_sheet_parameters'] = excel_sheet_parameters
+    excel_description['excel_sheet_structural_tower'] = excel_sheet_structural_tower
+    excel_description['excel_sheet_structural_blade'] = excel_sheet_structural_blade
+    excel_description['excel_sheet_discretization_blade'] = excel_sheet_discretization_blade
+    excel_description['excel_sheet_aero_blade'] = excel_sheet_aero_blade
+    excel_description['excel_sheet_airfoil_info'] = excel_sheet_airfoil_info
+    excel_description['excel_sheet_airfoil_chord'] = excel_sheet_airfoil_coord
+
+    wt = generate_from_excel_type03(op_params,
+                                   geom_params,
+                                   excel_description,
+                                   options)
+
+    return wt
diff --git a/lib/UVLM b/lib/UVLM
index eff039792..33a3fc86a 160000
--- a/lib/UVLM
+++ b/lib/UVLM
@@ -1,1 +1,1 @@
-Subproject commit eff039792d1e096564b89bb4f01975c798b0a7aa
+Subproject commit 33a3fc86a9b3835a1bc53c0272187229c4eb39f6
diff --git a/sharpy/aero/utils/uvlmlib.py b/sharpy/aero/utils/uvlmlib.py
index 9dd4b6d4f..06421014f 100644
--- a/sharpy/aero/utils/uvlmlib.py
+++ b/sharpy/aero/utils/uvlmlib.py
@@ -164,6 +164,7 @@ def vlm_solver(ts_info, options):
     flightconditions.uinf_direction = np.ctypeslib.as_ctypes(ts_info.u_ext[0][:, 0, 0]/flightconditions.uinf)
 
     p_rbm_vel_g = options['rbm_vel_g'].ctypes.data_as(ct.POINTER(ct.c_double))
+    p_centre_rot_g = options['centre_rot_g'].ctypes.data_as(ct.POINTER(ct.c_double))
 
     ts_info.generate_ctypes_pointers()
     run_VLM(ct.byref(vmopts),
@@ -177,7 +178,8 @@ def vlm_solver(ts_info, options):
             ts_info.ct_p_gamma,
             ts_info.ct_p_gamma_star,
             ts_info.ct_p_forces,
-            p_rbm_vel_g)
+            p_rbm_vel_g,
+            p_centre_rot_g)
     ts_info.remove_ctypes_pointers()
 
 
@@ -266,6 +268,7 @@ def uvlm_solver(i_iter, ts_info, struct_ts_info, options, convect_wake=True, dt=
     rbm_vel[0:3] = np.dot(struct_ts_info.cga(), rbm_vel[0:3])
     rbm_vel[3:6] = np.dot(struct_ts_info.cga(), rbm_vel[3:6])
     p_rbm_vel = rbm_vel.ctypes.data_as(ct.POINTER(ct.c_double))
+    p_centre_rot = options['centre_rot'].ctypes.data_as(ct.POINTER(ct.c_double))
 
     i = ct.c_uint(i_iter)
     ts_info.generate_ctypes_pointers()
@@ -281,6 +284,7 @@ def uvlm_solver(i_iter, ts_info, struct_ts_info, options, convect_wake=True, dt=
              ts_info.ct_p_zeta_star,
              ts_info.ct_p_zeta_dot,
              p_rbm_vel,
+             p_centre_rot,
              ts_info.ct_p_gamma,
              ts_info.ct_p_gamma_star,
              ts_info.ct_p_dist_to_orig,
diff --git a/sharpy/generators/floatingforces.py b/sharpy/generators/floatingforces.py
new file mode 100644
index 000000000..083b246b3
--- /dev/null
+++ b/sharpy/generators/floatingforces.py
@@ -0,0 +1,914 @@
+import numpy as np
+import h5py as h5
+import ctypes as ct
+import os
+from scipy import fft, ifft
+from scipy.interpolate import interp1d
+from control import forced_response, TransferFunction
+
+import sharpy.utils.cout_utils as cout
+import sharpy.utils.generator_interface as generator_interface
+import sharpy.utils.settings as settings
+import sharpy.utils.solver_interface as solver_interface
+from sharpy.utils.constants import deg2rad
+import sharpy.utils.h5utils as h5utils
+import sharpy.utils.algebra as algebra
+
+
+def compute_xf_zf(hf, vf, l, w, EA, cb):
+    """
+        Fairlead location (xf, zf) computation
+    """
+
+    root1, root2, ln1, ln2, lb = rename_terms(vf, hf, w, l)
+
+    # Define if there is part of the mooring line on the bed
+    if lb <= 0:
+        nobed = True
+    else:
+        nobed = False
+
+    # Compute the position of the fairlead
+    if nobed:
+        xf = hf/w*(ln1 - ln2) + hf*l/EA
+        zf = hf/w*(root1 - root2) + 1./EA*(vf*l-w*l**2/2)
+    else:
+        xf = lb + hf/w*ln1 + hf*l/EA
+        if not cb == 0.:
+            xf += cb*w/2/EA*(-lb**2 + (lb - hf/cb/w)*np.maximum((lb - hf/cb/w), 0))
+        zf = hf/w*(root1 - 1) + vf**2/2/EA/w
+    
+    return xf, zf
+
+
+def compute_jacobian(hf, vf, l, w, EA, cb):
+    """
+        Analytical computation of the Jacobian of equations
+        in function compute_xf_zf
+    """
+
+    root1, root2, ln1, ln2, lb = rename_terms(vf, hf, w, l)
+
+    # Compute their deivatives
+    der_root1_hf = 0.5*(1. + (vf/hf)**2)**(-0.5)*(2*vf/hf*(-vf/hf/hf))
+    der_root1_vf = 0.5*(1. + (vf/hf)**2)**(-0.5)*(2*vf/hf/hf)
+
+    der_root2_hf = 0.5*(1. + ((vf - w*l)/hf)**2)**(-0.5)*(2.*(vf - w*l)/hf*(-(vf - w*l)/hf/hf))
+    der_root2_vf = 0.5*(1. + ((vf - w*l)/hf)**2)**(-0.5)*(2.*(vf - w*l)/hf/hf)
+
+    der_ln1_hf = 1./(vf/hf + root1)*(vf/hf/hf + der_root1_hf)
+    der_ln1_vf = 1./(vf/hf + root1)*(1./hf + der_root1_vf)
+
+    der_ln2_hf = 1./((vf - w*l)/hf + root2)*(-(vf - w*l)/hf/hf + der_root2_hf)
+    der_ln2_vf = 1./((vf - w*l)/hf + root2)*(1./hf + der_root2_vf)
+
+    der_lb_hf = 0.
+    der_lb_vf = -1./w
+
+    # Define if there is part of the mooring line on the bed
+    if lb <= 0:
+        nobed = True
+    else:
+        nobed = False
+
+    # Compute the Jacobian
+    if nobed:
+        der_xf_hf = 1./w*(ln1 - ln2) + hf/w*(der_ln1_hf + der_ln2_hf) + l/EA
+        der_xf_vf = hf/w*(der_ln1_vf + der_ln2_vf)
+
+        der_zf_hf = 1./w*(root1 - root2) + hf/w*(der_root1_hf - der_root2_hf)
+        der_zf_vf = hf/w*(der_root1_vf - der_root2_vf) + 1./EA*l
+    else:
+        der_xf_hf = der_lb_hf + 1./w*ln1 + hf/w*der_ln1_hf + l/EA
+        if not cb == 0.:
+            arg1_max = l - vf/w - hf/cb/w
+            if arg1_max > 0.:
+                der_xf_hf += cb*w/2/EA*(2*(arg1_max)*(-1/cb/w))
+
+        der_xf_vf = der_lb_vf + hf/w*der_ln1_vf + cb*w/2/EA*(-2.*lb*der_lb_vf)
+        if not cb == 0.:
+            arg1_max = l - vf/w - hf/cb/w
+            if arg1_max > 0.:
+                der_xf_vf += cb*w/2/EA*(2.*(lb - hf/cb/w)*der_lb_vf)
+
+        der_zf_hf = 1/w*(root1 - 1) + hf/w*der_root1_hf
+        der_zf_vf = hf/w*der_root1_vf + vf/EA/w
+
+    J = np.array([[der_xf_hf, der_xf_vf],[der_zf_hf, der_zf_vf]])
+
+    return J
+
+
+def rename_terms(vf, hf, w, l):
+    """
+        Rename some terms for convenience
+    """
+    root1 = np.sqrt(1. + (vf/hf)**2)
+    root2 = np.sqrt(1. + ((vf - w*l)/hf)**2)
+    ln1 = np.log(vf/hf + root1)
+    ln2 = np.log((vf - w*l)/hf + root2)
+    lb = l - vf/w
+    return root1, root2, ln1, ln2, lb
+
+def quasisteady_mooring(xf, zf, l, w, EA, cb, hf0=None, vf0=None):
+    """
+        Computation of the forces generated by the mooring system
+        It performs a Newton-Raphson iteration based on the known equations
+        in compute_xf_zf function and the Jacobian
+    """
+
+    # Initialise guess for hf0 and vf0
+    if xf == 0:
+        lambda0 = 1e6
+    elif np.sqrt(xf**2 + zf**2) > l:
+        lambda0 = 0.2
+    else:
+        lambda0 = np.sqrt(3*((l**2 - zf**2)/xf**2 - 1))
+
+    if hf0 is None:
+        hf0 = np.abs(w*xf/2/lambda0)
+
+    if vf0 is None:
+        vf0 = w/2*(zf/np.tanh(lambda0) + l)
+
+    # Compute the solution through Newton-Raphson iteration
+    hf_est = hf0 + 0.
+    vf_est = vf0 + 0.
+    xf_est, zf_est = compute_xf_zf(hf_est, vf_est, l, w, EA, cb)
+    # print("initial: ", xf_est, zf_est)
+    tol = 1e-6
+    error = 2*tol
+    max_iter = 10000
+    it = 0
+    while ((error > tol) and (it < max_iter)):
+        J_est = compute_jacobian(hf_est, vf_est, l, w, EA, cb)
+        inv_J_est = np.linalg.inv(J_est)
+        hf_est += inv_J_est[0, 0]*(xf - xf_est) + inv_J_est[0, 1]*(zf - zf_est)
+        vf_est += inv_J_est[1, 0]*(xf - xf_est) + inv_J_est[1, 1]*(zf - zf_est)
+        # hf += (xf - xf_est)/J[0, 0] + (zf - zf_est)/J[1, 0]
+        # vf += (xf - xf_est)/J[0, 1] + (zf - zf_est)/J[1, 1]
+
+        xf_est, zf_est = compute_xf_zf(hf_est, vf_est, l, w, EA, cb)
+        error = np.maximum(np.abs(xf - xf_est), np.abs(zf - zf_est))
+        # print(error)
+        it += 1
+    if ((it == max_iter) and (error > tol)):
+        cout.cout_wrap(("Mooring system did not converge. error %f" % error), 4)
+        print("Mooring system did not converge. error %f" % error)
+
+    return hf_est, vf_est
+
+
+def wave_radiation_damping(K, qdot, it, dt):
+    """
+        This function computes the wave radiation damping assuming K constant
+    """
+    qdot_int = np.zeros((6,))
+    for idof in range(6):
+        qdot_int[idof] = np.trapz(np.arange(0, it + 1, 1)*dt, qdot[0:it, idof])
+
+    return np.dot(K, qdot_int)
+
+
+def change_of_to_sharpy(matrix_of):
+    """
+    Change between frame of reference of OpenFAST and the
+    usual one in SHARPy
+    """
+
+    sub_mat = np.array([[0., 0, 1],
+                        [0., -1, 0],
+                        [1., 0, 0]])
+    C_of_s = np.zeros((6,6))
+    C_of_s[0:3, 0:3] = sub_mat
+    C_of_s[3:6, 3:6] = sub_mat
+
+    matrix_sharpy = np.dot(C_of_s.T, np.dot(matrix_of, C_of_s))
+    return matrix_sharpy
+
+
+# def interp_1st_dim_matrix(A, vec, value):
+#
+#     # Make sure vec is ordered in strictly ascending order
+#     if (np.diff(vec) <= 0).any():
+#         cout.cout_wrap("ERROR: vec should be in strictly increasing order", 4)
+#     if not A.shape[0] == vec.shape[0]:
+#         cout.cout_wrap("ERROR: Incoherent vector and matrix size", 4)
+#
+#     # Compute the positions to interpolate
+#     if value <= vec[0]:
+#         return A[0, ...]
+#     elif ((value >= vec[-1]) or (value > vec[-2] and np.isinf(vec[-1]))):
+#         return A[-1, ...]
+#     else:
+#         i = 0
+#         while value > vec[i]:
+#             i += 1
+#         dist = vec[i] - vec[i - 1]
+#         rel_dist_to_im1 = (value - vec[i - 1])/dist
+#         rel_dist_to_i = (vec[i] - value)/dist
+#
+#     return A[i - 1, ...]*rel_dist_to_i + A[i, ...]*rel_dist_to_im1
+
+
+def rfval(num, den, z):
+    """
+        Evaluate a rational function given by the coefficients of the numerator (num) and
+        denominator (den) at z
+    """
+    return np.polyval(num, z)/np.polyval(den, z)
+
+
+def matrix_from_rf(dict_rf, w):
+    """
+    Create a matrix from the rational function approximation of each one of the elements
+    """
+    H = np.zeros((6, 6))
+    for i in range(6):
+        for j in range(6):
+            pos = "%d_%d" % (i, j)
+            H[i, j] = rfval(dict_rf[pos]['num'], dict_rf[pos]['den'], w)
+
+    return H
+
+
+def response_freq_dep_matrix(H, omega_H, q, it_, dt):
+    """
+    Compute the frequency response of a system with a transfer function depending on the frequency
+    F(t) = H(omega) * q(t)
+    """
+    it = it_ + 1
+    omega_fft = np.linspace(0, 1/(2*dt), it//2)[:it//2]
+    fourier_q = fft(q[:it, :], axis=0)
+    fourier_f = np.zeros_like(fourier_q)
+
+    ndof = q.shape[1]
+    f = np.zeros((ndof))
+
+    # Compute the constant component
+    if type(H) is np.ndarray:
+        # H_omega = interp_1st_dim_matrix(H, omega_H, omega_fft[0])
+        interp_H = interp1d(omega_H, H, axis=0)
+        H_omega = interp_H(omega_fft[0])
+    elif type(H) is tuple:
+        H_omega = matrix_from_rf(H, omega_fft[0])
+    else:
+        cout.cout_wrap(("ERROR: Not implemented response_freq_dep_matrix for type(H) %s" % type(H)), 4)
+    fourier_f[0, :] = np.dot(H_omega, fourier_q[0, :])
+
+    # Compute the rest of the terms
+    for iomega in range(1, omega_fft.shape[0]):
+        # Interpolate H at omega
+        if type(H) is np.ndarray:
+            H_omega = interp_H(omega_fft[iomega])
+        elif type(H) is dict:
+            H_omega = matrix_from_rf(H, omega_fft[iomega])
+        fourier_f[iomega, :] = np.dot(H_omega, fourier_q[iomega, :])
+        fourier_f[-iomega, :] = np.dot(H_omega, fourier_q[-iomega, :])
+
+    # Compute the inverse Fourier tranform
+    f[:] = np.real(ifft(fourier_f, axis=0)[it_, :])
+
+    # (T, yout, xout) = lsim(H, q[:it_ + 1, :], T, X0=X0)
+
+    return f
+
+
+def compute_equiv_hd_added_mass(f, q):
+    """
+        Compute the matrix H that satisfies f = Hq
+        H represents the added mass effects so it has to be
+        symmetric.
+        For the OC3 platfrom the following statements hold:
+            - z-y symmetry
+            - Non-diagonal non-zero terms: (1,5) and (2,4). Zero-indexed
+    """
+
+    if (q == 0).all():
+        return np.zeros((6,6))
+
+    q_mat = np.array([[q[0], 0,    0,    0,    0,    0],
+                      [0,    q[1], 0,    0,    q[5], 0],
+                      [0,    q[2], 0,    0,    0,    q[4]],
+                      [0,    0,    q[3], 0,    0,    0],
+                      [0,    0,    0,    q[4], 0,    q[2]],
+                      [0,    0,    0,    q[5], q[1], 0]])
+
+    hv = np.dot(np.linalg.inv(q_mat), f)
+
+    H = np.array([[hv[0], 0,     0,     0,     0,     0],
+                  [0,     hv[1], 0,     0,     0,     hv[4]],
+                  [0,     0,     hv[1], 0,     hv[5], 0],
+                  [0,     0,     0,     hv[2], 0,     0],
+                  [0,     0,     hv[5], 0,     hv[3], 0],
+                  [0,     hv[4], 0,     0,     0,     hv[3]]])
+
+    return H
+
+
+def jonswap_spectrum(Tp, Hs, w):
+    """
+    This function computes the one-sided spectrum of the JONSWAP wave data
+    [2] Jonkman, J. M. Dynamics modeling and loads analysis of an offshore floating wind turbine. 2007. NREL/TP-500-41958
+    """
+    nomega = w.shape[0]
+    spectrum = np.zeros((nomega))
+    for iomega in range(nomega):
+        # Compute the scaling factor
+        if w[iomega] <= 2*np.pi/Tp:
+            sigma = 0.07
+        else:
+            sigma = 0.09
+        # Compute the peak shape parameter
+        param = Tp/np.sqrt(Hs)
+        if param <= 3.6:
+            gamma = 5.
+        elif param > 5:
+            gamma = 1.
+        else:
+            gamma = np.exp(5.75 - 1.15*param)
+        # Compute one-sided spectrum
+        omega = w[iomega]
+        if omega == 0:
+            spectrum[iomega] = 0.
+        else:
+            param = omega*Tp/2/np.pi
+            spectrum[iomega] = (1./2/np.pi)*(5./16)*(Hs**2*Tp)*param**(-5)
+            spectrum[iomega] *= np.exp(-5./4*param**(-4))
+            spectrum[iomega] *= (1. - 0.287*np.log(gamma))
+            spectrum[iomega] *= gamma**np.exp(-0.5*((param - 1.)/sigma)**2)
+
+    return spectrum
+
+
+def noise_freq_1s(w):
+    """
+    Generates a frequency representation of a white noise
+    """
+    sigma = 1. #/np.sqrt(2)
+    nomega = w.shape[0]
+
+    wn = np.zeros((nomega, ), dtype=np.complex)
+    u1 = np.random.random(size=nomega) #+ 0j
+    u2 = np.random.random(size=nomega) #+ 0j
+    wn[0] = 0. + 0j
+    for iomega in range(1, nomega):
+        u1w = u1[iomega]
+        u2w = u2[iomega]
+        wn[iomega] = np.sqrt(-2.*np.log(u1w))*(np.cos(2*np.pi*u2w) +
+                             1j*np.sin(2*np.pi*u2w))
+        # wn[iomega] = 1. + 0j
+    return wn*sigma
+
+
+def time_wave_forces(Tp, Hs, dt, time, xi, w_xi):
+    """
+    Compute the time evolution of wave forces
+    """
+    # Compute time and frequency discretisations
+    ntime_steps = time.shape[0]
+
+    nomega = ntime_steps//2 + 1
+    w = np.zeros((nomega))
+    w_temp = np.fft.fftfreq(ntime_steps, d=dt)*2.*np.pi
+    w[:ntime_steps//2] = w_temp[:ntime_steps//2]
+    if ntime_steps%2 == 0:
+        w[-1] = -1*w_temp[ntime_steps//2]
+    else:
+        w[-1] = w_temp[ntime_steps//2]
+    nomega_2s = ntime_steps
+
+    # Compute the one-sided spectrums
+    noise_freq = noise_freq_1s(w)
+    jonswap_1s = jonswap_spectrum(Tp, Hs, w)*2.*np.pi
+    jonswap_freq = np.sqrt(2*ntime_steps/dt*jonswap_1s/2) + 0j
+    jonswap_freq[0] = np.sqrt(ntime_steps/dt*jonswap_1s[0]/2) + 0j # The DC values does not have the 2
+
+    # Compute the two-sided spectrum
+    force_freq_2s = np.zeros((nomega_2s, 6), dtype=np.complex)
+    for idim in range(6):
+        for iomega in range(nomega):
+            xi_interp = np.interp(w[iomega], w_xi, xi[:, idim])
+            force_freq_2s[iomega, idim] = (noise_freq[iomega]*
+                                           0.5*jonswap_freq[iomega]*
+                                           xi_interp)
+            if not iomega == 0:
+                if not ((iomega == nomega - 1) and (nomega_2s%2 == 0)):
+                    force_freq_2s[-iomega, idim] = (np.conj(noise_freq[iomega])*
+                                                    0.5*jonswap_freq[iomega]*
+                                                    np.conj(xi_interp))
+
+    # Compute the inverse Fourier transform
+    force_waves = ifft(force_freq_2s, axis=0)
+
+    return force_waves
+
+
+@generator_interface.generator
+class FloatingForces(generator_interface.BaseGenerator):
+    r"""
+    Floating forces generator
+
+    Generates the forces associated the floating support of offshore wind turbines.
+    Currently supports spar configurations.
+
+    The hydrostatic forces model includes buoyancy: an initial vertical force and the restoring forces
+    associated with heave, roll and pitch. See the implementation in [1].
+
+    The mooring model is the quasisteady implementation of Jonkman [2] .However, equation 2-37b is thought to be wrong (it is just a copy from eq 2-35b)
+    This was corrected according to the theory review of MAP++ [3].
+
+    The default values have been obtained from the OC3 platform report [1]
+
+    [1] Jonkman, J. Definition of the Floating System for Phase IV of OC3. 2010. NREL/TP-500-47535
+
+    [2] Jonkman, J. M. Dynamics modeling and loads analysis of an offshore floating wind turbine. 2007. NREL/TP-500-41958
+
+    [3] https://map-plus-plus.readthedocs.io/en/latest/theory.html (accessed on Octorber 14th, 2020)
+    """
+    generator_id = 'FloatingForces'
+
+    settings_types = dict()
+    settings_default = dict()
+    settings_description = dict()
+    settings_options = dict()
+
+    settings_types['n_time_steps'] = 'int'
+    settings_default['n_time_steps'] = None
+    settings_description['n_time_steps'] = 'Number of time steps'
+
+    settings_types['dt'] = 'float'
+    settings_default['dt'] = None
+    settings_description['dt'] = 'Time step'
+
+    settings_types['water_density'] = 'float'
+    settings_default['water_density'] = 1025 # kg/m3
+    settings_description['water_density'] = 'Water density'
+
+    settings_types['gravity'] = 'float'
+    settings_default['gravity'] = 9.81
+    settings_description['gravity'] = 'Gravity'
+
+    settings_types['gravity_dir'] = 'list(float)'
+    settings_default['gravity_dir'] = [1., 0., 0.]
+    settings_description['gravity_dir'] = 'Gravity direction'
+
+    settings_types['floating_file_name'] = 'str'
+    settings_default['floating_file_name'] = './oc3.floating.h5'
+    settings_description['floating_file_name'] = 'File containing the information about the floating dynamics'
+
+    settings_types['method_matrices_freq'] = 'str'
+    settings_default['method_matrices_freq'] = 'constant'
+    settings_description['method_matrices_freq'] = 'Method to compute frequency-dependent matrices'
+    settings_options['method_matrices_freq'] = ['constant', 'rational_function']
+
+    settings_types['matrices_freq'] = 'float'
+    settings_default['matrices_freq'] = 4.8 # Close to the upper limit defined in the oc3 report
+    settings_description['matrices_freq'] = 'Frequency [rad/s] to interpolate frequency-dependent matrices'
+
+    settings_types['steps_constant_matrices'] = 'int'
+    settings_default['steps_constant_matrices'] = 8
+    settings_description['steps_constant_matrices'] = 'Time steps to compute with constant matrices computed at ``matrices_freq``. Irrelevant in ``method_matrices_freq``=``constant``'
+
+    settings_types['added_mass_in_mass_matrix'] = 'bool'
+    settings_default['added_mass_in_mass_matrix'] = True
+    settings_description['added_mass_in_mass_matrix'] = 'Include the platform added mass in the mass matrix of the system'
+
+    settings_types['method_wave'] = 'str'
+    settings_default['method_wave'] = 'sin'
+    settings_description['method_wave'] = 'Method to compute wave forces'
+    settings_options['method_wave'] = ['sin', 'jonswap']
+
+    settings_types['wave_amplitude'] = 'float'
+    settings_default['wave_amplitude'] = 0.
+    settings_description['wave_amplitude'] = 'Wave amplitude. Only used in ``method_wave = sin``'
+
+    settings_types['wave_freq'] = 'float'
+    settings_default['wave_freq'] = 0.
+    settings_description['wave_freq'] = 'Wave circular frequency [rad/s]. Only used in ``method_wave = sin``'
+
+    settings_types['wave_Tp'] = 'float'
+    settings_default['wave_Tp'] = 0.
+    settings_description['wave_Tp'] = 'Wave peak spectral period [s]. Only used in ``method_wave = jonswap``'
+
+    settings_types['wave_Hs'] = 'float'
+    settings_default['wave_Hs'] = 0.
+    settings_description['wave_Hs'] = 'Significant wave height [m]. Only used in ``method_wave = jonswap``'
+
+    settings_types['wave_incidence'] = 'float'
+    settings_default['wave_incidence'] = 0.
+    settings_description['wave_incidence'] = 'Wave incidence in rad'
+
+    settings_types['write_output'] = 'bool'
+    settings_default['write_output'] = False
+    settings_description['write_output'] = 'Write forces to an output file'
+
+    settings_types['folder'] = 'str'
+    settings_default['folder'] = 'output'
+    settings_description['folder'] = 'Folder for the output files'
+
+    settings_types['log_filename'] = 'str'
+    settings_default['log_filename'] = 'log_floating_forces'
+    settings_description['log_filename'] = 'Log file name to write outputs'
+
+    setting_table = settings.SettingsTable()
+    __doc__ += setting_table.generate(settings_types, settings_default, settings_description)
+
+    def __init__(self):
+        self.in_dict = dict()
+
+        self.water_density = None
+        self.gravity = None
+        self.gravity_dir = None
+
+        self.floating_data = None
+
+        self.mooring_node = None
+        self.n_mooring_lines = None
+        self.anchor_pos = None
+        self.fairlead_pos_A = None
+        self.hf_prev = list() # Previous value of hf just for initialisation
+        self.vf_prev = list()
+
+        self.buoyancy_node = None
+        self.buoy_F0 = None
+        self.buoy_rest_mat = None
+
+        self.wave_forces_node = None
+
+        self.q = None
+        self.qdot = None
+        self.qdotdot = None
+
+        self.log_filename = None
+        self.added_mass_in_mass_matrix = None
+
+
+    def initialise(self, in_dict=None, data=None):
+        self.in_dict = in_dict
+        settings.to_custom_types(self.in_dict,
+                                 self.settings_types,
+                                 self.settings_default,
+                                 self.settings_options,
+                                 no_ctype=True)
+        self.settings = self.in_dict
+
+        self.water_density = self.settings['water_density']
+        self.gravity = self.settings['gravity']
+        self.gravity_dir = self.settings['gravity_dir']
+
+        # Platform dofs
+        self.q = np.zeros((self.settings['n_time_steps'] + 1, 6))
+        self.qdot = np.zeros_like(self.q)
+        self.qdotdot = np.zeros_like(self.q)
+
+        # Read the file with the floating information
+        fid = h5.File(self.settings['floating_file_name'], 'r')
+        self.floating_data = h5utils.load_h5_in_dict(fid)
+        fid.close()
+
+        # Mooringlines parameters
+        self.mooring_node = self.floating_data['mooring']['node']
+        self.n_mooring_lines = self.floating_data['mooring']['n_lines']
+        self.anchor_pos = np.zeros((self.n_mooring_lines, 3))
+        self.fairlead_pos_A = np.zeros((self.n_mooring_lines, 3))
+        self.hf_prev = [None]*self.n_mooring_lines
+        self.vf_prev = [None]*self.n_mooring_lines
+
+        theta = 2.*np.pi/self.n_mooring_lines
+        R = algebra.rotation3d_x(theta)
+        self.anchor_pos[0, 0] = -self.floating_data['mooring']['anchor_depth']
+        self.anchor_pos[0, 2] = self.floating_data['mooring']['anchor_radius']
+        self.fairlead_pos_A[0, 0] = -self.floating_data['mooring']['fairlead_depth']
+        self.fairlead_pos_A[0, 2] = self.floating_data['mooring']['fairlead_radius']
+        for imoor in range(1, self.n_mooring_lines):
+            self.anchor_pos[imoor, :] = np.dot(R, self.anchor_pos[imoor - 1, :])
+            self.fairlead_pos_A[imoor, :] = np.dot(R, self.fairlead_pos_A[imoor - 1, :])
+
+        # Hydrostatics
+        self.buoyancy_node = self.floating_data['hydrostatics']['node']
+        self.buoy_F0 = np.zeros((6,), dtype=float)
+        self.buoy_F0[0:3] = (self.floating_data['hydrostatics']['V0']*
+                              self.settings['water_density']*
+                              self.settings['gravity']*self.settings['gravity_dir'])
+        self.buoy_rest_mat = self.floating_data['hydrostatics']['buoyancy_restoring_matrix']
+
+        # hydrodynamics
+        if self.settings['method_matrices_freq'] == 'constant':
+            # self.hd_added_mass_const = interp_1st_dim_matrix(self.floating_data['hydrodynamics']['added_mass_matrix'],
+            #                             self.floating_data['hydrodynamics']['ab_freq_rads'],
+            #                             self.settings['matrices_freq'])
+
+            # self.hd_damping_const = interp_1st_dim_matrix(self.floating_data['hydrodynamics']['damping_matrix'],
+            #                             self.floating_data['hydrodynamics']['ab_freq_rads'],
+            #                             self.settings['matrices_freq'])
+
+            interp_am = interp1d(self.floating_data['hydrodynamics']['ab_freq_rads'],
+                                 self.floating_data['hydrodynamics']['added_mass_matrix'],
+                                 axis=0)
+            self.hd_added_mass_const = interp_am(self.settings['matrices_freq'])
+            interp_d = interp1d(self.floating_data['hydrodynamics']['ab_freq_rads'],
+                               self.floating_data['hydrodynamics']['damping_matrix'],
+                               axis=0)
+            self.hd_damping_const = interp_d(self.settings['matrices_freq'])
+
+
+        elif self.settings['method_matrices_freq'] == 'rational_function':
+            self.hd_added_mass_const = self.floating_data['hydrodynamics']['added_mass_matrix'][-1, :, :]
+            self.hd_damping_const = self.floating_data['hydrodynamics']['damping_matrix'][-1, :, :]
+
+        self.added_mass_in_mass_matrix = self.settings['added_mass_in_mass_matrix']
+        if self.added_mass_in_mass_matrix:
+            data.structure.add_lumped_mass_to_element(self.buoyancy_node,
+                                                      self.hd_added_mass_const)
+            data.structure.generate_fortran()
+
+        if self.settings['method_matrices_freq'] == 'rational_function':
+            ninput = 6
+            noutput = 6
+            hd_K_num = [None]*noutput
+            hd_K_den = [None]*noutput
+            for ioutput in range(noutput):
+                hd_K_num[ioutput] = [None]*ninput
+                hd_K_den[ioutput] = [None]*ninput
+                for iinput in range(ninput):
+                    pos = "%d_%d" % (ioutput, iinput)
+                    hd_K_num[ioutput][iinput] = self.floating_data['hydrodynamics']['K_rf'][pos]['num']
+                    hd_K_den[ioutput][iinput] = self.floating_data['hydrodynamics']['K_rf'][pos]['den']
+
+            self.hd_K = TransferFunction(hd_K_num, hd_K_den)
+            self.ab_freq_rads = self.floating_data['hydrodynamics']['ab_freq_rads']
+
+            self.x0_K = [None]*(self.settings['n_time_steps'] + 1)
+            self.x0_K[0] = 0.
+
+
+        # Wave forces
+        self.wave_forces_node = self.floating_data['wave_forces']['node']
+        if self.settings['method_wave'] == 'sin':
+            # xi_matrix2 = interp_1st_dim_matrix(self.floating_data['wave_forces']['xi'],
+            #                                 self.floating_data['wave_forces']['xi_freq_rads'],
+            #                                 self.settings['wave_freq'])
+            # xi = interp_1st_dim_matrix(xi_matrix2,
+            #                                 self.floating_data['wave_forces']['xi_beta_deg']*deg2rad,
+            #                                 self.settings['wave_incidence'])
+
+            interp_x1 = interp1d(self.floating_data['wave_forces']['xi_freq_rads'],
+                                 self.floating_data['wave_forces']['xi'],
+                                 axis=0,
+                                 bounds_error=False,
+                                 fill_value=(self.floating_data['wave_forces']['xi'][0, ...],
+                                             self.floating_data['wave_forces']['xi'][-1, ...]))
+            xi_matrix2 = interp_x1(self.settings['wave_freq'])
+            # xi_matrix2 = interp_x1(self.floating_data['wave_forces']['xi_freq_rads'][2:4])
+            interp_x2 = interp1d(self.floating_data['wave_forces']['xi_beta_deg']*deg2rad,
+                                 xi_matrix2,
+                                 axis=0)
+            xi = interp_x2(self.settings['wave_incidence'])
+            # xi = interp_x2(self.floating_data['wave_forces']['xi_beta_deg'][5]*deg2rad)
+            # print(xi[0, :])
+            # print(self.floating_data['wave_forces']['xi'][2, 5, :])
+            # print(xi[0, :] - self.floating_data['wave_forces']['xi'][2, 5, :])
+
+            phase = self.settings['wave_freq']*np.arange(self.settings['n_time_steps'] + 1)*self.settings['dt']
+            self.wave_forces_g = np.zeros((self.settings['n_time_steps'] + 1, 6))
+            for idim in range(6):
+                self.wave_forces_g[:, idim] = np.real(self.settings['wave_amplitude']*xi[idim]*(np.cos(phase) + 1j*np.sin(phase)))
+
+        elif self.settings['method_wave'] == 'jonswap':
+
+            interp_x1 = interp1d(self.floating_data['wave_forces']['xi_beta_deg']*deg2rad,
+                                 self.floating_data['wave_forces']['xi'],
+                                 axis=1)
+            xi_matrix = interp_x1(self.settings['wave_incidence'])
+
+            self.wave_forces_g = np.real(time_wave_forces(self.settings['wave_Tp'],
+                                                  self.settings['wave_Hs'],
+                                                  self.settings['dt'],
+                                                  np.arange(self.settings['n_time_steps'] + 1)*self.settings['dt'],
+                                                  xi_matrix,
+                                                  self.floating_data['wave_forces']['xi_freq_rads']))
+
+        # Log file
+        if not os.path.exists(self.settings['folder']):
+            os.makedirs(self.settings['folder'])
+        folder = self.settings['folder'] + '/' + data.settings['SHARPy']['case'] + '/floatingforces/'
+        if not os.path.exists(folder):
+            os.makedirs(folder)
+        self.log_filename = folder + self.settings['log_filename'] + '.h5'
+
+
+    def write_output(self, ts, k, mooring, mooring_yaw, hydrostatic,
+                     hydrodynamic_qdot, hydrodynamic_qdotdot, hd_correct_grav, waves):
+
+        output = dict()
+        output['ts'] = ts
+        output['k'] = k
+        output['q'] = self.q[ts, :]
+        output['qdot'] = self.qdot[ts, :]
+        output['qdotdot'] = self.qdotdot[ts, :]
+        output['mooring_forces'] = mooring
+        output['mooring_yaw'] = mooring_yaw
+        output['hydrostatic'] = hydrostatic
+        output['hydrodynamic_qdot'] = hydrodynamic_qdot
+        output['hydrodynamic_qdotdot'] = hydrodynamic_qdotdot
+        output['hydrodynamic_correct_grav'] = hd_correct_grav
+        output['waves'] = waves
+
+        fid = h5.File(self.log_filename, 'a')
+        group_name = "ts%d_k%d" % (ts, k)
+        if fid.__contains__(group_name):
+            del fid[group_name]
+        group = fid.create_group(group_name)
+        for key, value in output.items():
+            group.create_dataset(key, data=value)
+        fid.close()
+
+        debug_output = False
+        if debug_output:
+            print("q: ", self.q[ts, :])
+            print("qdot: ", self.qdot[ts, :])
+            print("qdotdot: ", self.qdotdot[ts, :])
+
+            print("mooring: ", mooring)
+            print("mooring_yaw: ", mooring_yaw)
+            print("hydrostatic: ", hydrostatic)
+            print("hydrodynamic_qdot: ", hydrodynamic_qdot)
+            print("hydrodynamic_qdotdot: ", hydrodynamic_qdotdot)
+            print("hydrodynamic_correct_grav: ", hd_correct_grav)
+            print("waves: ", waves)
+
+        return
+
+
+    def update_dof_vector(self, beam, struct_tstep, it, k):
+
+        if True:
+            cga = struct_tstep.cga()
+            self.q[it, 0:3] = (np.dot(cga, struct_tstep.pos[self.buoyancy_node, :]) +
+                      struct_tstep.for_pos[0:3])
+            self.q[it, 3:6] = algebra.quat2euler(struct_tstep.quat)
+
+            self.qdot[it, 0:3] = np.dot(cga, struct_tstep.for_vel[0:3])
+            self.qdot[it, 3:6] = np.dot(cga, struct_tstep.for_vel[3:6])
+
+            self.qdotdot[it, 0:3] = np.dot(cga, struct_tstep.for_acc[0:3])
+            self.qdotdot[it, 3:6] = np.dot(cga, struct_tstep.for_acc[3:6])
+
+        else:
+            self.q[it, :] = self.q[it-1, :]
+            self.qdot[it, :] = self.qdot[it-1, :]
+            self.qdotdot[it, :] = self.qdotdot[it-1, :]
+
+        return
+
+
+    def generate(self, params):
+        # Renaming for convenience
+        data = params['data']
+        struct_tstep = params['struct_tstep']
+        force_coeff = params['force_coeff']
+        k = params['fsi_substep']
+
+        # Update dof vector
+        self.update_dof_vector(data.structure, struct_tstep, data.ts, k)
+
+        # Mooring lines
+        mooring_forces = np.zeros((self.n_mooring_lines, 2))
+        cga = struct_tstep.cga()
+        ielem, inode_in_elem = data.structure.node_master_elem[self.mooring_node]
+        cab = algebra.crv2rotation(struct_tstep.psi[ielem, inode_in_elem])
+        cbg = np.dot(cab.T, cga.T)
+
+        moor_force_out = 0.
+        moor_mom_out = 0.
+
+        for imoor in range(self.n_mooring_lines):
+            fairlead_pos_G = (np.dot(cga, self.fairlead_pos_A[imoor, :]) +
+                              struct_tstep.for_pos[0:3])
+            fl_to_anchor_G = self.anchor_pos[imoor, :] - fairlead_pos_G
+            xf = np.sqrt(fl_to_anchor_G[1]**2 + fl_to_anchor_G[2]**2)
+            zf = np.abs(fl_to_anchor_G[0])
+            hf, vf = quasisteady_mooring(xf,
+                                 zf,
+                                 self.floating_data['mooring']['unstretched_length'],
+                                 self.floating_data['mooring']['apparent_weight'],
+                                 self.floating_data['mooring']['EA'],
+                                 self.floating_data['mooring']['seabed_drag_coef'],
+                                 hf0=self.hf_prev[imoor],
+                                 vf0=self.vf_prev[imoor])
+            mooring_forces[imoor, :] = np.array([hf, vf])
+            # Save the results to initialise the computation in the next time step
+            self.hf_prev[imoor] = hf + 0.
+            self.vf_prev[imoor] = vf + 0.
+
+            # Convert to the adequate reference system
+            horizontal_unit_vec = algebra.unit_vector(fl_to_anchor_G)
+            horizontal_unit_vec[0] = 0.
+            horizontal_unit_vec = algebra.unit_vector(horizontal_unit_vec)
+            force_fl = hf*horizontal_unit_vec + vf*np.array([-1., 0., 0.])
+
+            # Move the forces to the mooring node
+            force_cl = np.zeros((6,))
+            force_cl[0:3] = force_fl
+            mooring_node_pos_G = (np.dot(cga, struct_tstep.pos[self.mooring_node, :]) +
+                              struct_tstep.for_pos[0:3])
+            r_fairlead_G = fairlead_pos_G - mooring_node_pos_G
+            force_cl[3:6] = np.cross(r_fairlead_G, force_fl)
+
+            struct_tstep.runtime_generated_forces[self.mooring_node, 0:3] += np.dot(cbg, force_cl[0:3])
+            struct_tstep.runtime_generated_forces[self.mooring_node, 3:6] += np.dot(cbg, force_cl[3:6])
+
+        # Yaw moment generated by the mooring system
+        yaw = np.array([self.q[data.ts, 3], 0., 0.])
+        mooring_yaw = -self.floating_data['mooring']['yaw_spring_stif']*yaw
+        struct_tstep.runtime_generated_forces[self.mooring_node, 3:6] += np.dot(cbg,
+                                                                      mooring_yaw)
+
+        # Hydrostatic model
+        hs_f_g = self.buoy_F0 - np.dot(self.buoy_rest_mat, self.q[data.ts, :])
+
+        if not force_coeff == 0.:
+            hd_f_qdot_g = -np.dot(self.floating_data['hydrodynamics']['additional_damping'], self.qdot[data.ts, :])
+
+            if ((self.settings['method_matrices_freq'] == 'constant') or
+                 (data.ts < self.settings['steps_constant_matrices'])):
+                hd_f_qdot_g -= np.dot(self.hd_damping_const, self.qdot[data.ts, :])
+                hd_f_qdotdot_g = np.zeros((6))
+
+            elif self.settings['method_matrices_freq'] == 'rational_function':
+                # Damping
+                (T, yout, xout) = forced_response(self.hd_K,
+                                                  T=[0, self.settings['dt']],
+                                                  U=self.qdot[data.ts-1:data.ts+1, :].T,
+                                                  X0=self.x0_K[data.ts-1])
+                                                  # transpose=True)
+                self.x0_K[data.ts] = xout[:, 1]
+                hd_f_qdot_g -= yout[:, 1]
+                hd_f_qdotdot_g = np.zeros((6))
+
+            else:
+                cout.cout_wrap(("ERROR: Unknown method_matrices_freq %s" % self.settings['method_matrices_freq']), 4)
+
+            # Correct gravity forces if needed
+            if self.added_mass_in_mass_matrix:
+                # Correct unreal gravity forces from added mass
+                gravity_b = np.zeros((6,),)
+                gravity_b[0:3] = np.dot(cbg, -self.settings['gravity_dir'])*self.settings['gravity']
+                hd_correct_grav = -np.dot(self.hd_added_mass_const, gravity_b)
+                struct_tstep.runtime_generated_forces[self.buoyancy_node, :] += hd_correct_grav
+            else:
+                hd_correct_grav = np.zeros((6))
+
+        else:
+            hd_f_qdot_g = np.zeros((6))
+            hd_f_qdotdot_g = np.zeros((6))
+            hd_correct_grav = np.zeros((6))
+
+        ielem, inode_in_elem = data.structure.node_master_elem[self.buoyancy_node]
+        cab = algebra.crv2rotation(struct_tstep.psi[ielem, inode_in_elem])
+        cbg = np.dot(cab.T, cga.T)
+        struct_tstep.runtime_generated_forces[self.buoyancy_node, 0:3] += np.dot(cbg, hs_f_g[0:3] + force_coeff*(hd_f_qdot_g[0:3] + hd_f_qdotdot_g[0:3]))
+        struct_tstep.runtime_generated_forces[self.buoyancy_node, 3:6] += np.dot(cbg, hs_f_g[3:6] + force_coeff*(hd_f_qdot_g[3:6] + hd_f_qdotdot_g[3:6]))
+
+        # Nonlinear drag coefficeint
+        for inode in range(self.floating_data['hydrodynamics']['CD_first_node'], self.floating_data['hydrodynamics']['CD_last_node'] + 1):
+
+            ielem, inode_in_elem = data.structure.node_master_elem[inode]
+            cab = algebra.crv2rotation(struct_tstep.psi[ielem, inode_in_elem])
+            cbg = np.dot(cab.T, cga.T)
+
+            if inode == 0:
+                delta_x = 0.5*np.linalg.norm(struct_tstep.pos[1, :] - struct_tstep.pos[0, :])
+            elif inode == data.structure.num_node - 1:
+                delta_x = 0.5*np.linalg.norm(struct_tstep.pos[inode, :] - struct_tstep.pos[inode - 1, :])
+            else:
+                delta_x = 0.5*np.linalg.norm(struct_tstep.pos[inode + 1, :] - struct_tstep.pos[inode - 1, :])
+
+            vel_a = (struct_tstep.for_vel[0:3] +
+                     np.cross(struct_tstep.for_vel[3:6], struct_tstep.pos[inode, :]) +
+                     struct_tstep.pos_dot[inode, :])
+
+            # Remove velocity along the x axis
+            vel_b = np.dot(cab.T, vel_a)
+            vel_b[0] = 0.
+            vel_g = np.dot(cbg.T, vel_b)
+
+            drag_force = (-0.5*self.water_density*np.linalg.norm(vel_g)*vel_g*delta_x*
+                          self.floating_data['hydrodynamics']['spar_diameter']*
+                          self.floating_data['hydrodynamics']['CD'])
+
+            struct_tstep.runtime_generated_forces[inode, 0:3] += np.dot(cbg, drag_force)
+
+        # Wave loading
+        ielem, inode_in_elem = data.structure.node_master_elem[self.wave_forces_node]
+        cab = algebra.crv2rotation(struct_tstep.psi[ielem, inode_in_elem])
+        cbg = np.dot(cab.T, cga.T)
+
+        struct_tstep.runtime_generated_forces[self.wave_forces_node, 0:3] += np.dot(cbg, force_coeff*self.wave_forces_g[data.ts, 0:3])
+        struct_tstep.runtime_generated_forces[self.wave_forces_node, 3:6] += np.dot(cbg, force_coeff*self.wave_forces_g[data.ts, 3:6])
+
+        # Write output
+        if self.settings['write_output']:
+            self.write_output(data.ts, k, mooring_forces, mooring_yaw, hs_f_g,
+                     hd_f_qdot_g, hd_f_qdotdot_g, hd_correct_grav, self.wave_forces_g[data.ts, :])
diff --git a/sharpy/postproc/savedata.py b/sharpy/postproc/savedata.py
index feb1f89a5..7cb9fab7c 100644
--- a/sharpy/postproc/savedata.py
+++ b/sharpy/postproc/savedata.py
@@ -6,11 +6,9 @@
 from sharpy.utils.solver_interface import solver, BaseSolver
 import sharpy.utils.settings as settings
 import sharpy.utils.h5utils as h5utils
+from sharpy.presharpy.presharpy import PreSharpy
 
 
-# Define basic numerical types
-# BasicNumTypes=(float,float32,float64,int,int32,int64,complex)
-
 @solver
 class SaveData(BaseSolver):
     """
@@ -55,6 +53,10 @@ class SaveData(BaseSolver):
     settings_description['save_linear_uvlm'] = 'Save linear UVLM state space system. Use with caution when dealing with ' \
                                                'large systems.'
 
+    settings_types['save_wake'] = 'bool'
+    settings_default['save_wake'] = True
+    settings_description['save_wake'] = 'Save aero wake classes.'
+
     settings_types['save_rom'] = 'bool'
     settings_default['save_rom'] = False
     settings_description['save_rom'] = 'Saves the ROM matrices and the reduced order model'
@@ -93,7 +95,6 @@ class SaveData(BaseSolver):
                                        settings_options=settings_options)
 
     def __init__(self):
-        import sharpy
 
         self.settings = None
         self.data = None
@@ -106,7 +107,7 @@ def __init__(self):
 
         ### specify which classes are saved as hdf5 group
         # see initialise and add_as_grp
-        self.ClassesToSave = (sharpy.presharpy.presharpy.PreSharpy,)
+        self.ClassesToSave = (PreSharpy,)
 
     def initialise(self, data, custom_settings=None, caller=None):
         self.data = data
@@ -164,6 +165,12 @@ def initialise(self, data, custom_settings=None, caller=None):
             if self.settings['save_aero']:
                 self.ClassesToSave += (sharpy.aero.models.aerogrid.Aerogrid,
                                        sharpy.utils.datastructures.AeroTimeStepInfo,)
+                if not self.settings['save_wake']:
+                    self.settings['skip_attr'].append('zeta_star')
+                    self.settings['skip_attr'].append('u_ext_star')
+                    self.settings['skip_attr'].append('gamma_star')
+                    self.settings['skip_attr'].append('dist_to_orig')
+                    self.settings['skip_attr'].append('wake_conv_vel')
 
             if self.settings['save_struct']:
                 self.ClassesToSave += (
diff --git a/sharpy/solvers/_basestructural.py b/sharpy/solvers/_basestructural.py
index 5534f66f2..e1bf1cbf1 100644
--- a/sharpy/solvers/_basestructural.py
+++ b/sharpy/solvers/_basestructural.py
@@ -50,6 +50,10 @@ class _BaseStructural(BaseSolver):
     settings_default['gravity'] = 9.81
     settings_description['gravity'] = 'Gravitational acceleration'
 
+    settings_types['gravity_dir'] = 'list(float)'
+    settings_default['gravity_dir'] = [0., 0., 1.]
+    settings_description['gravity_dir'] = 'Direction in G where gravity applies'
+
     settings_types['relaxation_factor'] = 'float'
     settings_default['relaxation_factor'] = 0.3
 
diff --git a/sharpy/solvers/beamloader.py b/sharpy/solvers/beamloader.py
index e13339022..95b92e66f 100644
--- a/sharpy/solvers/beamloader.py
+++ b/sharpy/solvers/beamloader.py
@@ -56,6 +56,10 @@ class BeamLoader(BaseSolver):
     settings_default['orientation'] = [1., 0, 0, 0]
     settings_description['orientation'] = 'Initial attitude of the structure given as the quaternion that parametrises the rotation from G to A frames of reference.'
 
+    settings_types['for_pos'] = 'list(float)'
+    settings_default['for_pos'] = [0., 0, 0]
+    settings_description['for_pos'] = 'Initial position of the A FoR.'
+
     settings_table = settings.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
diff --git a/sharpy/solvers/dynamiccoupled.py b/sharpy/solvers/dynamiccoupled.py
index 21e227771..204428d78 100644
--- a/sharpy/solvers/dynamiccoupled.py
+++ b/sharpy/solvers/dynamiccoupled.py
@@ -493,12 +493,26 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None):
                     structural_kstep, aero_kstep = self.process_controller_output(
                         state)
 
+            k = 0
+            # compute unsteady contribution
+            force_coeff = 0.0
+            unsteady_contribution = False
+            if self.settings['include_unsteady_force_contribution']:
+                if self.data.ts > self.settings['steps_without_unsteady_force']:
+                    unsteady_contribution = True
+                    if 0 < self.settings['pseudosteps_ramp_unsteady_force']:
+                        force_coeff = k/self.settings['pseudosteps_ramp_unsteady_force']
+                    else:
+                        force_coeff = 1.
             # Add external forces
             if self.with_runtime_generators:
+                structural_kstep.runtime_generated_forces.fill(0.)
                 params = dict()
                 params['data'] = self.data
                 params['struct_tstep'] = structural_kstep
                 params['aero_tstep'] = aero_kstep
+                params['force_coeff'] = force_coeff
+                params['fsi_substep'] = -1
                 for id, runtime_generator in self.runtime_generators.items():
                     runtime_generator.generate(params)
 
@@ -510,12 +524,12 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None):
             controlled_structural_kstep = structural_kstep.copy()
             controlled_aero_kstep = aero_kstep.copy()
 
-            k = 0
             for k in range(self.settings['fsi_substeps'] + 1):
                 if (k == self.settings['fsi_substeps'] and
                         self.settings['fsi_substeps']):
-                    print_res = 0 if self.res_dqdt == 0. else np.log10(self.res_dqdt)
-                    cout.cout_wrap(("The FSI solver did not converge!!! residual: %f" % print_res))
+                    print_res = 0 if self.res == 0. else np.log10(self.res)
+                    print_res_dqdt = 0 if self.res_dqdt == 0. else np.log10(self.res_dqdt)
+                    cout.cout_wrap(("The FSI solver did not converge!!! residuals: %f %f" % (print_res, print_res_dqdt)))
                     self.aero_solver.update_custom_grid(
                         structural_kstep,
                         aero_kstep)
@@ -538,6 +552,19 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None):
                         else:
                             force_coeff = 1.
 
+                previous_runtime_generated_forces = structural_kstep.runtime_generated_forces.astype(dtype=ct.c_double, order='F', copy=True)
+                # Add external forces
+                if self.with_runtime_generators:
+                    structural_kstep.runtime_generated_forces.fill(0.)
+                    params = dict()
+                    params['data'] = self.data
+                    params['struct_tstep'] = structural_kstep
+                    params['aero_tstep'] = aero_kstep
+                    params['force_coeff'] = force_coeff
+                    params['fsi_substep'] = k
+                    for id, runtime_generator in self.runtime_generators.items():
+                        runtime_generator.generate(params)
+
                 # run the solver
                 ini_time_aero = time.perf_counter()
                 self.data = self.aero_solver.run(aero_kstep,
@@ -548,10 +575,13 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None):
 
                 previous_kstep = structural_kstep.copy()
                 structural_kstep = controlled_structural_kstep.copy()
+                structural_kstep.runtime_generated_forces = previous_kstep.runtime_generated_forces.astype(dtype=ct.c_double, order='F', copy=True)
+                previous_kstep.runtime_generated_forces = previous_runtime_generated_forces.astype(dtype=ct.c_double, order='F', copy=True)
 
                 # move the aerodynamic surface according the the structural one
                 self.aero_solver.update_custom_grid(structural_kstep,
                                                     aero_kstep)
+
                 self.map_forces(aero_kstep,
                                 structural_kstep,
                                 force_coeff)
@@ -593,7 +623,10 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None):
                 # check convergence
                 if self.convergence(k,
                                     structural_kstep,
-                                    previous_kstep) or self.settings['aero_solver'].lower() == 'noaero':
+                                    previous_kstep,
+                                    self.settings['structural_solver'].lower(),
+                                    self.settings['aero_solver'].lower(),
+                                    self.with_runtime_generators):
                     # move the aerodynamic surface according to the structural one
                     self.aero_solver.update_custom_grid(
                         structural_kstep,
@@ -643,7 +676,8 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None):
             finish_event.set()
             self.logger.info('Time loop - Complete')
 
-    def convergence(self, k, tstep, previous_tstep):
+    def convergence(self, k, tstep, previous_tstep,
+                    struct_solver, aero_solver, with_runtime_generators):
         r"""
         Check convergence in the FSI loop.
 
@@ -668,8 +702,16 @@ def convergence(self, k, tstep, previous_tstep):
             self.base_dqdt = np.linalg.norm(tstep.dqdt.copy())
             if self.base_dqdt == 0:
                 self.base_dqdt = 1.
+            if with_runtime_generators:
+                self.base_res_forces = np.linalg.norm(tstep.runtime_generated_forces)
+                if self.base_res_forces == 0:
+                    self.base_res_forces = 1.
             return False
 
+        # Check the special case of no aero and no runtime generators
+        if aero_solver.lower() == "noaero" and not with_runtime_generators:
+            return True
+
         # relative residuals
         self.res = (np.linalg.norm(tstep.q-
                                    previous_tstep.q)/
@@ -678,6 +720,13 @@ def convergence(self, k, tstep, previous_tstep):
                                         previous_tstep.dqdt)/
                          self.base_dqdt)
 
+        if with_runtime_generators:
+            res_forces = (np.linalg.norm(tstep.runtime_generated_forces -
+                                        previous_tstep.runtime_generated_forces)/
+                         self.base_res_forces)
+        else:
+            res_forces = 0.
+
         # we don't want this to converge before introducing the gamma_dot forces!
         if self.settings['include_unsteady_force_contribution']:
             if k < self.settings['pseudosteps_ramp_unsteady_force'] \
@@ -685,12 +734,9 @@ def convergence(self, k, tstep, previous_tstep):
                 return False
 
         # convergence
-        if k > self.settings['minimum_steps'] - 1:
-            if self.res < self.settings['fsi_tolerance']:
-                if self.res_dqdt < self.settings['fsi_tolerance']:
-                    return True
+        return k > self.settings['minimum_steps'] - 1 and \
+            all(x < self.settings['fsi_tolerance'] for x in (self.res, self.res_dqdt, res_forces))
 
-        return False
 
     def map_forces(self, aero_kstep, structural_kstep, unsteady_forces_coeff=1.0):
         # set all forces to 0
@@ -731,18 +777,17 @@ def map_forces(self, aero_kstep, structural_kstep, unsteady_forces_coeff=1.0):
             #                                                      dynamic_struct_forces)
 
         # prescribed forces + aero forces
+        structural_kstep.steady_applied_forces = (
+            (struct_forces + self.data.structure.ini_info.steady_applied_forces).
+            astype(dtype=ct.c_double, order='F', copy=True))
         try:
-            structural_kstep.steady_applied_forces = (
-                (struct_forces + self.data.structure.ini_info.steady_applied_forces).
-                astype(dtype=ct.c_double, order='F', copy=True))
             structural_kstep.unsteady_applied_forces = (
-                (dynamic_struct_forces + self.data.structure.dynamic_input[max(self.data.ts - 1, 0)]['dynamic_forces']).
+                (dynamic_struct_forces + self.data.structure.dynamic_input[max(self.data.ts - 1, 0)]['dynamic_forces'] +
+                 structural_kstep.runtime_generated_forces).
                 astype(dtype=ct.c_double, order='F', copy=True))
         except KeyError:
-            structural_kstep.steady_applied_forces = (
-                (struct_forces + self.data.structure.ini_info.steady_applied_forces).
-                astype(dtype=ct.c_double, order='F', copy=True))
-            structural_kstep.unsteady_applied_forces = dynamic_struct_forces
+            structural_kstep.unsteady_applied_forces = (dynamic_struct_forces +
+                                                        structural_kstep.runtime_generated_forces)
 
     def relaxation_factor(self, k):
         initial = self.settings['relaxation_factor']
@@ -795,10 +840,12 @@ def interpolate_timesteps(step0, step1, out_step, coeff):
 
 
 def relax(beam, timestep, previous_timestep, coeff):
-    timestep.steady_applied_forces[:] = ((1.0 - coeff)*timestep.steady_applied_forces +
+    timestep.steady_applied_forces = ((1.0 - coeff)*timestep.steady_applied_forces +
             coeff*previous_timestep.steady_applied_forces)
-    timestep.unsteady_applied_forces[:] = ((1.0 - coeff)*timestep.unsteady_applied_forces +
+    timestep.unsteady_applied_forces = ((1.0 - coeff)*timestep.unsteady_applied_forces +
             coeff*previous_timestep.unsteady_applied_forces)
+    timestep.runtime_generated_forces = ((1.0 - coeff)*timestep.runtime_generated_forces +
+            coeff*previous_timestep.runtime_generated_forces)
 
 
 def normalise_quaternion(tstep):
diff --git a/sharpy/solvers/nonlineardynamicmultibody.py b/sharpy/solvers/nonlineardynamicmultibody.py
index 772208a33..a93abef41 100644
--- a/sharpy/solvers/nonlineardynamicmultibody.py
+++ b/sharpy/solvers/nonlineardynamicmultibody.py
@@ -84,6 +84,12 @@ def initialise(self, data, custom_settings=None):
         self.data.structure.add_unsteady_information(
             self.data.structure.dyn_dict, self.settings['num_steps'])
 
+        # Initialise time integrator
+        self.time_integrator = solver_interface.initialise_solver(
+            self.settings['time_integrator'])
+        self.time_integrator.initialise(
+            self.data, self.settings['time_integrator_settings'])
+
         # Define the number of equations
         self.lc_list = lagrangeconstraints.initialize_constraints(self.data.structure.ini_mb_dict)
         self.num_LM_eq = lagrangeconstraints.define_num_LM_eq(self.lc_list)
diff --git a/sharpy/solvers/staticcoupled.py b/sharpy/solvers/staticcoupled.py
index 422e9e422..6a1d429c7 100644
--- a/sharpy/solvers/staticcoupled.py
+++ b/sharpy/solvers/staticcoupled.py
@@ -9,6 +9,7 @@
 import sharpy.utils.settings as settings
 import sharpy.utils.algebra as algebra
 import sharpy.utils.correct_forces as cf
+import sharpy.utils.generator_interface as gen_interface
 
 @solver
 class StaticCoupled(BaseSolver):
@@ -65,6 +66,12 @@ class StaticCoupled(BaseSolver):
     settings_description['correct_forces_method'] = 'Function used to correct aerodynamic forces. Check :py:mod:`sharpy.utils.correct_forces`'
     settings_options['correct_forces_method'] = ['efficiency', 'polars']
 
+    settings_types['runtime_generators'] = 'dict'
+    settings_default['runtime_generators'] = dict()
+    settings_description['runtime_generators'] = 'The dictionary keys are the runtime generators to be used. ' \
+                                                 'The dictionary values are dictionaries with the settings ' \
+                                                 'needed by each generator.'
+
     settings_table = settings.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
 
@@ -82,6 +89,9 @@ def __init__(self):
         self.correct_forces = False
         self.correct_forces_function = None
 
+        self.runtime_generators = dict()
+        self.with_runtime_generators = False
+
     def initialise(self, data, input_dict=None):
         self.data = data
         if input_dict is None:
@@ -113,6 +123,15 @@ def initialise(self, data, input_dict=None):
             self.correct_forces = True
             self.correct_forces_function = cf.dict_of_corrections[self.settings['correct_forces_method']]
 
+        # initialise runtime generators
+        self.runtime_generators = dict()
+        if self.settings['runtime_generators']:
+            self.with_runtime_generators = True
+            for id, param in self.settings['runtime_generators'].items():
+                gen = gen_interface.generator_from_string(id)
+                self.runtime_generators[id] = gen()
+                self.runtime_generators[id].initialise(param, data=self.data)
+
     def increase_ts(self):
         self.data.ts += 1
         self.structural_solver.next_step()
@@ -169,6 +188,20 @@ def run(self):
                                         struct_forces,
                                         rho=self.aero_solver.settings['rho'])
 
+                # Add external forces
+                if self.with_runtime_generators:
+                    self.data.structure.timestep_info[self.data.ts].runtime_generated_forces.fill(0.)
+                    params = dict()
+                    params['data'] = self.data
+                    params['struct_tstep'] = self.data.structure.timestep_info[self.data.ts]
+                    params['aero_tstep'] = self.data.aero.timestep_info[self.data.ts]
+                    params['force_coeff'] = 0.
+                    params['fsi_substep'] = -i_iter
+                    for id, runtime_generator in self.runtime_generators.items():
+                        runtime_generator.generate(params)
+
+                    struct_forces += self.data.structure.timestep_info[self.data.ts].runtime_generated_forces
+
                 if not self.settings['relaxation_factor'] == 0.:
                     if i_iter == 0:
                         self.previous_force = struct_forces.copy()
diff --git a/sharpy/solvers/staticuvlm.py b/sharpy/solvers/staticuvlm.py
index bddcc34a2..f924525bb 100644
--- a/sharpy/solvers/staticuvlm.py
+++ b/sharpy/solvers/staticuvlm.py
@@ -106,6 +106,10 @@ class StaticUvlm(BaseSolver):
     settings_default['rbm_vel_g'] = [0., 0., 0., 0., 0., 0.]
     settings_description['rbm_vel_g'] = 'Rigid body velocity in G FoR'
 
+    settings_types['centre_rot_g'] = 'list(float)'
+    settings_default['centre_rot_g'] = [0., 0., 0.]
+    settings_description['centre_rot_g'] = 'Centre of rotation in G FoR around which ``rbm_vel_g`` is applied'
+
     settings_table = settings.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
diff --git a/sharpy/solvers/stepuvlm.py b/sharpy/solvers/stepuvlm.py
index 3d00d7442..566cc8325 100644
--- a/sharpy/solvers/stepuvlm.py
+++ b/sharpy/solvers/stepuvlm.py
@@ -125,6 +125,10 @@ class StepUvlm(BaseSolver):
     settings_default['yaw_slerp'] = 0
     settings_description['yaw_slerp'] = 'Yaw angle in radians to be used when interp_metod == 4'
 
+    settings_types['centre_rot'] = 'list(float)'
+    settings_default['centre_rot'] = [0., 0., 0.]
+    settings_description['centre_rot'] = 'Coordinates of the centre of rotation to perform slerp interpolation or cylindrical coordinates'
+
     settings_types['quasi_steady'] = 'bool'
     settings_default['quasi_steady'] = False
     settings_description['quasi_steady'] = 'Use quasi-steady approximation in UVLM'
diff --git a/sharpy/structure/models/beam.py b/sharpy/structure/models/beam.py
index 3c893166f..614531529 100644
--- a/sharpy/structure/models/beam.py
+++ b/sharpy/structure/models/beam.py
@@ -41,6 +41,8 @@ def __init__(self):
         self.lumped_mass_nodes = None
         self.lumped_mass_inertia = None
         self.lumped_mass_position = None
+        self.lumped_mass_mat = None
+        self.lumped_mass_mat_nodes = None
         self.n_lumped_mass = 0
 
         self.steady_app_forces = None
@@ -176,11 +178,13 @@ def generate(self, in_data, settings):
         self.ini_info.steady_applied_forces = self.steady_app_forces.astype(dtype=ct.c_double, order='F')
         # rigid body rotations
         self.ini_info.quat = self.settings['orientation'].astype(dtype=ct.c_double, order='F')
+        self.ini_info.for_pos[0:3] = self.settings['for_pos'].astype(dtype=ct.c_double, order='F')
 
         self.timestep_info.append(self.ini_info.copy())
         self.timestep_info[-1].steady_applied_forces = self.steady_app_forces.astype(dtype=ct.c_double, order='F')
 
         # lumped masses
+        self.n_lumped_mass = 0
         try:
             self.lumped_mass = in_data['lumped_mass'].copy()
         except KeyError:
@@ -189,9 +193,19 @@ def generate(self, in_data, settings):
             self.lumped_mass_nodes = in_data['lumped_mass_nodes'].copy()
             self.lumped_mass_inertia = in_data['lumped_mass_inertia'].copy()
             self.lumped_mass_position = in_data['lumped_mass_position'].copy()
-            self.n_lumped_mass, _ = self.lumped_mass_position.shape
+            self.n_lumped_mass += self.lumped_mass_position.shape[0]
+
+        # lumped masses given as matrices
+        try:
+            self.lumped_mass_mat = in_data['lumped_mass_mat'].copy()
+        except KeyError:
+            self.lumped_mass_mat = None
+        else:
+            self.lumped_mass_mat_nodes = in_data['lumped_mass_mat_nodes'].copy()
+            self.n_lumped_mass += self.lumped_mass_mat.shape[0]
+
         # lumped masses to element mass
-        if self.lumped_mass is not None:
+        if self.lumped_mass is not None or self.lumped_mass_mat is not None:
             self.lump_masses()
 
         # self.generate_dof_arrays()
@@ -272,31 +286,52 @@ def generate_dof_arrays(self):
 
         self.num_dof = ct.c_int((vcounter + 1)*6)
 
-    def lump_masses(self):
-        for i_lumped in range(self.n_lumped_mass):
-            r = self.lumped_mass_position[i_lumped, :]
-            m = self.lumped_mass[i_lumped]
-            j = self.lumped_mass_inertia[i_lumped, :, :]
+    def generate_mass_matrix(self, mass, position, inertia):
 
-            i_lumped_node = self.lumped_mass_nodes[i_lumped]
-            i_lumped_master_elem, i_lumped_master_node_local = self.node_master_elem[i_lumped_node]
+        inertia_tensor = np.zeros((6, 6))
+        pos_skew = algebra.skew(position)
+        inertia_tensor[0:3, 0:3] = mass*np.eye(3)
+        inertia_tensor[0:3, 3:6] = -mass*pos_skew
+        inertia_tensor[3:6, 0:3] = mass*pos_skew
+        inertia_tensor[3:6, 3:6] = inertia + mass*(np.dot(pos_skew.T, pos_skew))
+
+        return inertia_tensor
 
-            # cba = algebra.crv2rot(self.elements[i_lumped_master_elem].psi_def[i_lumped_master_node_local, :]).T
+    def lump_masses(self):
+        if self.lumped_mass is not None:
+            for i_lumped in range(self.lumped_mass.shape[0]):
+                r = self.lumped_mass_position[i_lumped, :]
+                m = self.lumped_mass[i_lumped]
+                j = self.lumped_mass_inertia[i_lumped, :, :]
+
+                inertia_tensor = self.generate_mass_matrix(m, r, j)
+                i_lumped_node = self.lumped_mass_nodes[i_lumped]
+                self.add_lumped_mass_to_element(i_lumped_node,
+                                                   inertia_tensor)
+        if self.lumped_mass_mat is not None:
+            for i_lumped in range(self.lumped_mass_mat.shape[0]):
+                inertia_tensor = self.lumped_mass_mat[i_lumped, :, :]
+                i_lumped_node = self.lumped_mass_mat_nodes[i_lumped]
+                self.add_lumped_mass_to_element(i_lumped_node,
+                                                   inertia_tensor)
+
+        return
+
+    def add_lumped_mass_to_element(self, i_lumped_node, inertia_tensor, replace=False):
 
-            inertia_tensor = np.zeros((6, 6))
-            r_skew = algebra.skew(r)
-            inertia_tensor[0:3, 0:3] = m*np.eye(3)
-            inertia_tensor[0:3, 3:6] = -m*r_skew
-            inertia_tensor[3:6, 0:3] = m*r_skew
-            inertia_tensor[3:6, 3:6] = j + m*(np.dot(r_skew.T, r_skew))
+            i_lumped_master_elem, i_lumped_master_node_local = self.node_master_elem[i_lumped_node]
 
             if self.elements[i_lumped_master_elem].rbmass is None:
                 # allocate memory
                 self.elements[i_lumped_master_elem].rbmass = np.zeros((
                     self.elements[i_lumped_master_elem].max_nodes_elem, 6, 6))
 
-            self.elements[i_lumped_master_elem].rbmass[i_lumped_master_node_local, :, :] += (
-                inertia_tensor) # += necessary in case multiple masses defined per node
+            if replace:
+                self.elements[i_lumped_master_elem].rbmass[i_lumped_master_node_local, :, :] = (
+                    inertia_tensor) # += necessary in case multiple masses defined per node
+            else:
+                self.elements[i_lumped_master_elem].rbmass[i_lumped_master_node_local, :, :] = (
+                    inertia_tensor) # += necessary in case multiple masses defined per node
 
     # def generate_master_structure(self):
     #     self.master = np.zeros((self.num_elem, self.num_node_elem, 2), dtype=int) - 1
diff --git a/sharpy/structure/utils/xbeamlib.py b/sharpy/structure/utils/xbeamlib.py
index 39a690f28..936685069 100644
--- a/sharpy/structure/utils/xbeamlib.py
+++ b/sharpy/structure/utils/xbeamlib.py
@@ -90,8 +90,7 @@ def cbeam3_solv_nlnstatic(beam, settings, ts):
     xbopts.MinDelta = ct.c_double(settings['min_delta'])
     xbopts.gravity_on = ct.c_bool(settings['gravity_on'])
     xbopts.gravity = ct.c_double(settings['gravity'])
-    gravity_vector = np.array([0.0, 0.0, 1.0])
-    gravity_vector = np.dot(beam.timestep_info[ts].cag(), gravity_vector)
+    gravity_vector = np.dot(beam.timestep_info[ts].cag(), settings['gravity_dir'])
     xbopts.gravity_dir_x = ct.c_double(gravity_vector[0])
     xbopts.gravity_dir_y = ct.c_double(gravity_vector[1])
     xbopts.gravity_dir_z = ct.c_double(gravity_vector[2])
@@ -164,13 +163,11 @@ def cbeam3_solv_nlndyn(beam, settings):
     f_cbeam3_solv_nlndyn = xbeamlib.cbeam3_solv_nlndyn_python
     f_cbeam3_solv_nlndyn.restype = None
 
-
     n_elem = ct.c_int(beam.num_elem)
     n_nodes = ct.c_int(beam.num_node)
     n_mass = ct.c_int(beam.n_mass)
     n_stiff = ct.c_int(beam.n_stiff)
 
-
     dt = settings['dt']
     n_tsteps = settings['num_steps']
     time = np.zeros((n_tsteps,), dtype=ct.c_double, order='F')
@@ -276,10 +273,10 @@ def cbeam3_step_nlndyn(beam, settings, ts, tstep=None, dt=None):
     xbopts.NewmarkDamp = ct.c_double(settings['newmark_damp'])
     xbopts.gravity_on = ct.c_bool(settings['gravity_on'])
     xbopts.gravity = ct.c_double(settings['gravity'])
-    xbopts.gravity_dir_x = ct.c_double(tstep.gravity_vector_inertial[0])
-    xbopts.gravity_dir_y = ct.c_double(tstep.gravity_vector_inertial[1])
-    xbopts.gravity_dir_z = ct.c_double(tstep.gravity_vector_inertial[2])
-    xbopts.relaxation_factor = settings['relaxation_factor']
+    xbopts.gravity_dir_x = ct.c_double(settings['gravity_dir'][0])
+    xbopts.gravity_dir_y = ct.c_double(settings['gravity_dir'][1])
+    xbopts.gravity_dir_z = ct.c_double(settings['gravity_dir'][2])
+    xbopts.relaxation_factor = ct.c_double(settings['relaxation_factor'])
 
     # here we only need to set the flags at True, all the forces are follower
     xbopts.FollowerForce = ct.c_bool(True)
@@ -370,9 +367,9 @@ def xbeam_solv_couplednlndyn(beam, settings):
     xbopts.NewmarkDamp = ct.c_double(settings['newmark_damp'])
     xbopts.gravity_on = ct.c_bool(settings['gravity_on'])
     xbopts.gravity = ct.c_double(settings['gravity'])
-    xbopts.gravity_dir_x = ct.c_double(beam.ini_info.gravity_vector_inertial[0])
-    xbopts.gravity_dir_y = ct.c_double(beam.ini_info.gravity_vector_inertial[1])
-    xbopts.gravity_dir_z = ct.c_double(beam.ini_info.gravity_vector_inertial[2])
+    xbopts.gravity_dir_x = ct.c_double(settings['gravity_dir'][0])
+    xbopts.gravity_dir_y = ct.c_double(settings['gravity_dir'][1])
+    xbopts.gravity_dir_z = ct.c_double(settings['gravity_dir'][2])
     xbopts.relaxation_factor = ct.c_double(settings['relaxation_factor'])
 
     pos_def_history = np.zeros((n_tsteps.value, beam.num_node, 3), order='F', dtype=ct.c_double)
@@ -478,9 +475,9 @@ def xbeam_step_couplednlndyn(beam, settings, ts, tstep=None, dt=None):
     xbopts.gravity_on = ct.c_bool(settings['gravity_on'])
     xbopts.gravity = ct.c_double(settings['gravity'])
     xbopts.balancing = ct.c_bool(settings['balancing'])
-    xbopts.gravity_dir_x = ct.c_double(tstep.gravity_vector_inertial[0])
-    xbopts.gravity_dir_y = ct.c_double(tstep.gravity_vector_inertial[1])
-    xbopts.gravity_dir_z = ct.c_double(tstep.gravity_vector_inertial[2])
+    xbopts.gravity_dir_x = ct.c_double(settings['gravity_dir'][0])
+    xbopts.gravity_dir_y = ct.c_double(settings['gravity_dir'][1])
+    xbopts.gravity_dir_z = ct.c_double(settings['gravity_dir'][2])
     xbopts.relaxation_factor = ct.c_double(settings['relaxation_factor'])
 
     if dt is None:
@@ -552,9 +549,9 @@ def xbeam_init_couplednlndyn(beam, settings, ts, dt=None):
     xbopts.NewmarkDamp = ct.c_double(settings['newmark_damp'])
     xbopts.gravity_on = ct.c_bool(settings['gravity_on'])
     xbopts.gravity = ct.c_double(settings['gravity'])
-    xbopts.gravity_dir_x = ct.c_double(beam.timestep_info[ts].gravity_vector_inertial[0])
-    xbopts.gravity_dir_y = ct.c_double(beam.timestep_info[ts].gravity_vector_inertial[1])
-    xbopts.gravity_dir_z = ct.c_double(beam.timestep_info[ts].gravity_vector_inertial[2])
+    xbopts.gravity_dir_x = ct.c_double(settings['gravity_dir'][0])
+    xbopts.gravity_dir_y = ct.c_double(settings['gravity_dir'][1])
+    xbopts.gravity_dir_z = ct.c_double(settings['gravity_dir'][2])
 
     in_dt = ct.c_double(dt)
 
@@ -886,9 +883,9 @@ def cbeam3_asbly_dynamic(beam, tstep, settings):
     xbopts.NewmarkDamp = ct.c_double(settings['newmark_damp'])
     xbopts.gravity_on = ct.c_bool(settings['gravity_on'])
     xbopts.gravity = ct.c_double(settings['gravity'])
-    xbopts.gravity_dir_x = ct.c_double(tstep.gravity_vector_inertial[0])
-    xbopts.gravity_dir_y = ct.c_double(tstep.gravity_vector_inertial[1])
-    xbopts.gravity_dir_z = ct.c_double(tstep.gravity_vector_inertial[2])
+    xbopts.gravity_dir_x = ct.c_double(settings['gravity_dir'][0])
+    xbopts.gravity_dir_y = ct.c_double(settings['gravity_dir'][1])
+    xbopts.gravity_dir_z = ct.c_double(settings['gravity_dir'][2])
 
     # Initialize matrices
     Mglobal = np.zeros((num_dof, num_dof), dtype=ct.c_double, order='F')
@@ -990,9 +987,9 @@ def xbeam3_asbly_dynamic(beam, tstep, settings):
     xbopts.NewmarkDamp = ct.c_double(settings['newmark_damp'])
     xbopts.gravity_on = ct.c_bool(settings['gravity_on'])
     xbopts.gravity = ct.c_double(settings['gravity'])
-    xbopts.gravity_dir_x = ct.c_double(tstep.gravity_vector_inertial[0])
-    xbopts.gravity_dir_y = ct.c_double(tstep.gravity_vector_inertial[1])
-    xbopts.gravity_dir_z = ct.c_double(tstep.gravity_vector_inertial[2])
+    xbopts.gravity_dir_x = ct.c_double(settings['gravity_dir'][0])
+    xbopts.gravity_dir_y = ct.c_double(settings['gravity_dir'][1])
+    xbopts.gravity_dir_z = ct.c_double(settings['gravity_dir'][2])
 
     # Initialize matrices
     Mtotal = np.zeros((num_dof+10, num_dof+10), dtype=ct.c_double, order='F')
@@ -1137,9 +1134,10 @@ def cbeam3_asbly_static(beam, tstep, settings, iLoadStep):
     # xbopts.NewmarkDamp = settings['newmark_damp']
     xbopts.gravity_on = ct.c_bool(settings['gravity_on'])
     xbopts.gravity = ct.c_double(settings['gravity'])
-    xbopts.gravity_dir_x = ct.c_double(tstep.gravity_vector_inertial[0])
-    xbopts.gravity_dir_y = ct.c_double(tstep.gravity_vector_inertial[1])
-    xbopts.gravity_dir_z = ct.c_double(tstep.gravity_vector_inertial[2])
+    gravity_vector = np.dot(beam.timestep_info[ts].cag(), settings['gravity_dir'])
+    xbopts.gravity_dir_x = ct.c_double(gravity_vector[0])
+    xbopts.gravity_dir_y = ct.c_double(gravity_vector[1])
+    xbopts.gravity_dir_z = ct.c_double(gravity_vector[2])
 
     # Initialize matrices
     # Mglobal = np.zeros((num_dof, num_dof), dtype=ct.c_double, order='F')
@@ -1204,10 +1202,10 @@ def xbeam_step_coupledrigid(beam, settings, ts, tstep=None, dt=None):
     xbopts.gravity_on = ct.c_bool(settings['gravity_on'])
     xbopts.gravity = ct.c_double(settings['gravity'])
     xbopts.balancing = ct.c_bool(settings['balancing'])
-    xbopts.gravity_dir_x = ct.c_double(tstep.gravity_vector_inertial[0])
-    xbopts.gravity_dir_y = ct.c_double(tstep.gravity_vector_inertial[1])
-    xbopts.gravity_dir_z = ct.c_double(tstep.gravity_vector_inertial[2])
-    xbopts.relaxation_factor = settings['relaxation_factor']
+    xbopts.gravity_dir_x = ct.c_double(settings['gravity_dir'][0])
+    xbopts.gravity_dir_y = ct.c_double(settings['gravity_dir'][1])
+    xbopts.gravity_dir_z = ct.c_double(settings['gravity_dir'][2])
+    xbopts.relaxation_factor = ct.c_double(settings['relaxation_factor'])
 
     if dt is None:
         in_dt = ct.c_double(settings['dt'])
diff --git a/sharpy/utils/datastructures.py b/sharpy/utils/datastructures.py
index 346313a3a..d07ee9776 100644
--- a/sharpy/utils/datastructures.py
+++ b/sharpy/utils/datastructures.py
@@ -503,13 +503,12 @@ class StructTimeStepInfo(object):
         for_vel (np.ndarray): ``A`` frame of reference velocity. Expressed in A FoR
         for_acc (np.ndarray): ``A`` frame of reference acceleration. Expressed in A FoR
 
-        gravity_vector_inertial (np.ndarray): Unit vector parallel to the direction of gravity expressed in the ``G``
-          FoR
-        gravity_vector_body (np.ndarray): Unit vector parallel to the direction of gravity expressed in the ``A`` FoR
         steady_applied_forces (np.ndarray): Forces applied to the structure not associated to time derivatives
           ``[num_nodes x 6]``. Expressed in B FoR
         unsteady_applied_forces (np.ndarray): Forces applied to the structure associated to time derivatives
           ``[num_node x 6]``. Expressed in B FoR
+        runtime_generated_forces (np.ndarray): Forces generated at runtime through runtime generators
+          ``[num_node x 6]``. Expressed in B FoR
         gravity_forces (np.ndarray): Gravity forces at nodes ``[num_node x 6]``. Expressed in A FoR
 
         total_gravity_forces (np.ndarray): Total gravity forces on the structure ``[6]``. Expressed in A FoR
@@ -555,11 +554,9 @@ def __init__(self, num_node, num_elem, num_node_elem=3, num_dof=None, num_bodies
         self.for_vel = np.zeros((6,), dtype=ct.c_double, order='F')
         self.for_acc = np.zeros((6,), dtype=ct.c_double, order='F')
 
-        self.gravity_vector_inertial = np.array([0.0, 0.0, 1.0], dtype=ct.c_double, order='F')
-        self.gravity_vector_body = np.array([0.0, 0.0, 1.0], dtype=ct.c_double, order='F')
-
         self.steady_applied_forces = np.zeros((self.num_node, 6), dtype=ct.c_double, order='F')
         self.unsteady_applied_forces = np.zeros((self.num_node, 6), dtype=ct.c_double, order='F')
+        self.runtime_generated_forces = np.zeros((self.num_node, 6), dtype=ct.c_double, order='F')
         self.gravity_forces = np.zeros((self.num_node, 6), dtype=ct.c_double, order='F')
         self.total_gravity_forces = np.zeros((6,), dtype=ct.c_double, order='F')
         self.total_forces = np.zeros((6,), dtype=ct.c_double, order='F')
@@ -612,11 +609,9 @@ def copy(self):
         copied.for_vel = self.for_vel.astype(dtype=ct.c_double, order='F', copy=True)
         copied.for_acc = self.for_acc.astype(dtype=ct.c_double, order='F', copy=True)
 
-        copied.gravity_vector_inertial = self.gravity_vector_inertial.astype(dtype=ct.c_double, order='F', copy=True)
-        copied.gravity_vector_body = self.gravity_vector_body.astype(dtype=ct.c_double, order='F', copy=True)
-
         copied.steady_applied_forces = self.steady_applied_forces.astype(dtype=ct.c_double, order='F', copy=True)
         copied.unsteady_applied_forces = self.unsteady_applied_forces.astype(dtype=ct.c_double, order='F', copy=True)
+        copied.runtime_generated_forces = self.runtime_generated_forces.astype(dtype=ct.c_double, order='F', copy=True)
         copied.gravity_forces = self.gravity_forces.astype(dtype=ct.c_double, order='F', copy=True)
         copied.total_gravity_forces = self.total_gravity_forces.astype(dtype=ct.c_double, order='F', copy=True)
         copied.total_forces = self.total_forces.astype(dtype=ct.c_double, order='F', copy=True)
@@ -715,11 +710,9 @@ def get_body(self, beam, num_dof_ibody, ibody):
         ibody_StructTimeStepInfo.psi_dot = self.psi_dot[ibody_elems,:,:].astype(dtype=ct.c_double, order='F', copy=True)
         ibody_StructTimeStepInfo.psi_ddot = self.psi_ddot[ibody_elems,:,:].astype(dtype=ct.c_double, order='F', copy=True)
 
-        ibody_StructTimeStepInfo.gravity_vector_inertial = self.gravity_vector_inertial.astype(dtype=ct.c_double, order='F', copy=True)
-        ibody_StructTimeStepInfo.gravity_vector_body = self.gravity_vector_body.astype(dtype=ct.c_double, order='F', copy=True)
-
         ibody_StructTimeStepInfo.steady_applied_forces = self.steady_applied_forces[ibody_nodes,:].astype(dtype=ct.c_double, order='F', copy=True)
         ibody_StructTimeStepInfo.unsteady_applied_forces = self.unsteady_applied_forces[ibody_nodes,:].astype(dtype=ct.c_double, order='F', copy=True)
+        ibody_StructTimeStepInfo.runtime_generated_forces = self.runtime_generated_forces[ibody_nodes,:].astype(dtype=ct.c_double, order='F', copy=True)
         ibody_StructTimeStepInfo.gravity_forces = self.gravity_forces[ibody_nodes,:].astype(dtype=ct.c_double, order='F', copy=True)
         ibody_StructTimeStepInfo.total_gravity_forces = self.total_gravity_forces.astype(dtype=ct.c_double, order='F', copy=True)
 
diff --git a/sharpy/utils/generate_cases.py b/sharpy/utils/generate_cases.py
index 042a91361..9d36be63f 100644
--- a/sharpy/utils/generate_cases.py
+++ b/sharpy/utils/generate_cases.py
@@ -167,6 +167,11 @@ def read_column_sheet_type01(excel_file_name, excel_sheet, column_name):
     excel_db = pd.read_excel(xls, sheet_name=excel_sheet)
     num_elem = excel_db.index.stop - 2
 
+    try:
+        excel_db[column_name]
+    except KeyError:
+        return None
+
     if excel_db[column_name][1] == 'one_int':
         var = excel_db[column_name][2]
     elif excel_db[column_name][1] == 'one_float':
@@ -346,6 +351,9 @@ def __init__(self):
         self.lumped_mass = None
         self.lumped_mass_inertia = None
         self.lumped_mass_position = None
+        self.lumped_mass_mat = None
+        self.lumped_mass_mat_nodes = None
+
 
     def copy(self):
         """
@@ -378,11 +386,15 @@ def copy(self):
             copied.lumped_mass = self.lumped_mass.astype(dtype=float, copy=True)
             copied.lumped_mass_inertia = self.lumped_mass_inertia.astype(dtype=float, copy=True)
             copied.lumped_mass_position = self.lumped_mass_position.astype(dtype=float, copy=True)
+        if isinstance(self.lumped_mass_mat_nodes, np.ndarray):
+            copied.lumped_mass_mat_nodes = self.lumped_mass_mat_nodes.astype(dtype=int, copy=True)
+            copied.lumped_mass_mat = self.lumped_mass_mat.astype(dtype=float, copy=True)
 
         return copied
 
     def set_to_zero(self, num_node_elem, num_node, num_elem,
-                    num_mass_db=None, num_stiffness_db=None, num_lumped_mass=0):
+                    num_mass_db=None, num_stiffness_db=None, num_lumped_mass=0,
+                    num_lumped_mass_mat=0):
         """
         set_to_zero
 
@@ -395,6 +407,7 @@ def set_to_zero(self, num_node_elem, num_node, num_elem,
             num_mass_db (int): number of different mass matrices in the case
             num_stiffness_db (int): number of different stiffness matrices in the case
             num_lumped_mass (int): number of lumped masses in the case
+            num_lumped_mass_mat (int): number of lumped masses given as matrices
         """
 
         if num_mass_db is None:
@@ -425,6 +438,10 @@ def set_to_zero(self, num_node_elem, num_node, num_elem,
                                                 dtype=float)
             self.lumped_mass_position = np.zeros((num_lumped_mass, 3),
                                                  dtype=float)
+        if not num_lumped_mass_mat == 0:
+            self.lumped_mass_mat_nodes = np.zeros((num_lumped_mass_mat,), dtype=int)
+            self.lumped_mass_mat = np.zeros((num_lumped_mass_mat, 6, 6), dtype=float)
+
 
     def generate_full_structure(self,
                                 num_node_elem,
@@ -444,7 +461,9 @@ def generate_full_structure(self,
                                 lumped_mass_nodes=None,
                                 lumped_mass=None,
                                 lumped_mass_inertia=None,
-                                lumped_mass_position=None):
+                                lumped_mass_position=None,
+                                lumped_mass_mat_nodes=None,
+                                lumped_mass_mat=None):
         """
         generate_full_structure
 
@@ -469,6 +488,8 @@ def generate_full_structure(self,
             lumped_mass (np.array): value of the lumped masses
             lumped_mass_inertia (np.array): inertia of the lumped masses
             lumped_mass_position (np.array): position of the lumped masses
+            lumped_mass_mat_nodes (np.array): nodes with lumped masses given by matrices
+            lumped_mass_mat (np.array): value of the lumped masses given by matrices
         """
 
         self.num_node_elem = num_node_elem
@@ -491,6 +512,10 @@ def generate_full_structure(self,
             self.lumped_mass = lumped_mass
             self.lumped_mass_inertia = lumped_mass_inertia
             self.lumped_mass_position = lumped_mass_position
+        if isinstance(self.lumped_mass_mat_nodes, np.ndarray):
+            self.lumped_mass_mat_nodes = lumped_mass_mat_nodes
+            self.lumped_mass_mat = lumped_mass_mat
+
 
     def generate_1to1_from_vectors(self,
                                     num_node_elem,
@@ -501,7 +526,8 @@ def generate_1to1_from_vectors(self,
                                     mass_db,
                                     frame_of_reference_delta,
                                     vec_node_structural_twist,
-                                    num_lumped_mass=0):
+                                    num_lumped_mass=0,
+                                    num_lumped_mass_mat=0):
 
         self.num_node_elem = num_node_elem
         self.num_node = num_node
@@ -514,6 +540,7 @@ def generate_1to1_from_vectors(self,
         self.mass_db = mass_db
         self.create_frame_of_reference_delta(y_BFoR=frame_of_reference_delta)
         self.structural_twist = from_node_list_to_elem_matrix(vec_node_structural_twist, self.connectivities)
+        self.boundary_conditions = np.zeros((self.num_node,), dtype=int)
         self.beam_number = np.zeros((self.num_elem,), dtype=int)
         self.body_number = np.zeros((self.num_elem,), dtype=int)
         self.app_forces = np.zeros((self.num_node, 6), dtype=float)
@@ -522,6 +549,10 @@ def generate_1to1_from_vectors(self,
             self.lumped_mass = np.zeros((num_lumped_mass,), dtype=float)
             self.lumped_mass_inertia = np.zeros((num_lumped_mass, 3, 3), dtype=float)
             self.lumped_mass_position = np.zeros((num_lumped_mass, 3), dtype=float)
+        if not num_lumped_mass_mat == 0:
+            self.lumped_mass_mat_nodes = np.zeros((num_lumped_mass_mat,), dtype=int)
+            self.lumped_mass_mat = np.zeros((num_lumped_mass_mat, 6, 6), dtype=float)
+
 
     def create_frame_of_reference_delta(self, y_BFoR='y_AFoR'):
         """
@@ -686,7 +717,8 @@ def generate_uniform_sym_beam(self,
                                   EI,
                                   num_node_elem=3,
                                   y_BFoR='y_AFoR',
-                                  num_lumped_mass=0):
+                                  num_lumped_mass=0,
+                                  num_lumped_mass_mat=0):
         """
         generate_uniform_sym_beam
 
@@ -718,7 +750,8 @@ def generate_uniform_sym_beam(self,
                                    EI,
                                    num_node_elem,
                                    y_BFoR,
-                                   num_lumped_mass)
+                                   num_lumped_mass,
+                                   num_lumped_mass_mat)
 
     def generate_uniform_beam(self,
                               node_pos,
@@ -735,7 +768,8 @@ def generate_uniform_beam(self,
                               EIz,
                               num_node_elem=3,
                               y_BFoR='y_AFoR',
-                              num_lumped_mass=0):
+                              num_lumped_mass=0,
+                              num_lumped_mass_mat=0):
         """
         generate_uniform_beam
 
@@ -757,12 +791,14 @@ def generate_uniform_beam(self,
             num_node_elem (int): number of nodes per element
             y_BFoR (str): orientation of the yB axis
             num_lumped_mass (int): number of lumped masses
+            num_lumped_mass_mat (int): number of lumped masses given as matrices
         """
         self.num_node = len(node_pos)
         self.num_node_elem = num_node_elem
         self.compute_basic_num_elem()
 
-        self.set_to_zero(self.num_node_elem, self.num_node, self.num_elem, 1, 1, num_lumped_mass)
+        self.set_to_zero(self.num_node_elem, self.num_node, self.num_elem, 1, 1,
+                         num_lumped_mass, num_lumped_mass_mat)
         self.coordinates = node_pos
         self.create_simple_connectivities()
         # self.create_mass_db_from_vector(np.ones((num_elem,),)*mass_per_unit_length,
@@ -790,6 +826,7 @@ def generate_uniform_beam(self,
         self.create_frame_of_reference_delta(y_BFoR)
         self.boundary_conditions = np.zeros((self.num_node), dtype=int)
 
+
     def assembly_structures(self, *args):
         """
         assembly_structures
@@ -811,6 +848,7 @@ def assembly_structures(self, *args):
         total_num_mass = self.mass_db.shape[0]
 
         for structure_to_add in args:
+
             self.coordinates = np.concatenate((self.coordinates, structure_to_add.coordinates), axis=0)
             self.connectivities = np.concatenate((self.connectivities, structure_to_add.connectivities + total_num_node), axis=0)
             assert self.num_node_elem == structure_to_add.num_node_elem, "num_node_elem does NOT match"
@@ -835,6 +873,12 @@ def assembly_structures(self, *args):
                 self.lumped_mass = structure_to_add.lumped_mass
                 self.lumped_mass_inertia = structure_to_add.lumped_mass_inertia
                 self.lumped_mass_position = structure_to_add.lumped_mass_position
+            if isinstance(self.lumped_mass_mat_nodes, np.ndarray) and isinstance(structure_to_add.lumped_mass_mat_nodes, np.ndarray):
+                self.lumped_mass_mat_nodes = np.concatenate((self.lumped_mass_mat_nodes, structure_to_add.lumped_mass_mat_nodes + total_num_node), axis=0)
+                self.lumped_mass_mat = np.concatenate((self.lumped_mass_mat, structure_to_add.lumped_mass_mat), axis=0)
+            elif isinstance(structure_to_add.lumped_mass_mat_nodes, np.ndarray):
+                self.lumped_mass_mat_nodes = structure_to_add.lumped_mass_mat_nodes + total_num_node
+                self.lumped_mass_mat = structure_to_add.lumped_mass_mat
 
             total_num_stiff += structure_to_add.stiffness_db.shape[0]
             total_num_mass += structure_to_add.mass_db.shape[0]
@@ -919,6 +963,9 @@ def generate_fem_file(self, route, case_name):
                 h5file.create_dataset('lumped_mass', data=self.lumped_mass)
                 h5file.create_dataset('lumped_mass_inertia', data=self.lumped_mass_inertia)
                 h5file.create_dataset('lumped_mass_position', data=self.lumped_mass_position)
+            if isinstance(self.lumped_mass_mat_nodes, np.ndarray):
+                h5file.create_dataset('lumped_mass_mat_nodes', data=self.lumped_mass_mat_nodes)
+                h5file.create_dataset('lumped_mass_mat', data=self.lumped_mass_mat)
 
 
 ######################################################################
@@ -1442,7 +1489,7 @@ def assembly(self, *args):
         self.StructuralInformation.assembly_structures(*list_of_SI)
         self.AerodynamicInformation.assembly_aerodynamics(*list_of_AI)
 
-    def remove_duplicated_points(self, tol):
+    def remove_duplicated_points(self, tol, skip=[]):
         """
         remove_duplicated_points
 
@@ -1450,6 +1497,7 @@ def remove_duplicated_points(self, tol):
 
         Args:
             tol (float): tolerance. Maximum distance between nodes to be merged
+            skip (list): nodes to keep (do not remove)
 
         Notes:
             This function will not work if an element or an aerdoynamic surface is completely eliminated
@@ -1475,9 +1523,10 @@ def find_connectivities_index(connectivities, iprev_node):
         # Fill the first three columns
         for inode in range(self.StructuralInformation.num_node):
             for iprev_node in range(inode):
-                if (np.linalg.norm(
+                if ((np.linalg.norm(
                     self.StructuralInformation.coordinates[inode, :] -
-                    self.StructuralInformation.coordinates[iprev_node, :]) < tol):
+                    self.StructuralInformation.coordinates[iprev_node, :]) < tol) and (
+                       inode not in skip)):
                     cout.cout_wrap(("WARNING: Replacing node %d by node %d" % (inode, iprev_node)), 3)
                     replace_matrix[inode, 0] = iprev_node
                     replace_matrix[inode, 1], replace_matrix[inode, 2] = (
@@ -1643,35 +1692,32 @@ def set_default_values(self):
 
         # TODO: I am sure this can be done in a better way
         for solver in aux_names:
-            self.solvers[solver] = {}
-            if not solver == 'PreSharpy':
-                try:
-                    aux_solver = solver_interface.solver_from_string(solver)
-                except:
-                    aux_solver = generator_interface.generator_from_string(solver)
-                    # TODO: remove this try/except when generators are rewriten as solvers with class attributes instead of instance attributes
-                    aux_solver.__init__(aux_solver)
+            if solver == 'PreSharpy':
+                solver_name = 'SHARPy'
+            else:
+                solver_name = solver
+            self.solvers[solver_name] = {}
+            try:
+                aux_solver = solver_interface.solver_from_string(solver)
+            except:
+                aux_solver = generator_interface.generator_from_string(solver)
+                # TODO: remove this try/except when generators are rewriten as solvers with class attributes instead of instance attributes
+                aux_solver.__init__(aux_solver)
 
-                if aux_solver.settings_default == {}:
-                    aux_solver.__init__(aux_solver)
+            if aux_solver.settings_default == {}:
+                aux_solver.__init__(aux_solver)
 
-                for key, value in aux_solver.settings_default.items():
-                    self.solvers[solver][key] = deepcopy(value)
+            for key, value in aux_solver.settings_default.items():
+                self.solvers[solver_name][key] = deepcopy(value)
 
         # # MAIN
-        self.solvers['SHARPy'] = {'flow': '',
-                                  'case': 'default_case_name',
-                                  'route': '',
-                                  'write_screen': 'on',
-                                  'write_log': 'off',
-                                  'log_folder': './output',
-                                  'log_file': 'log'}
-
-        self.solvers['SaveData'] = {'folder': './output',
-                                   'save_aero': True,
-                                   'save_struct': True,
-                                   # 'skip_attr': dict(),
-                                   'compress_float': False}
+        # self.solvers['SHARPy'] = {'flow': '',
+        #                           'case': 'default_case_name',
+        #                           'route': '',
+        #                           'write_screen': 'on',
+        #                           'write_log': 'off',
+        #                           'log_folder': './output',
+        #                           'log_file': 'log'}
 
     def check(self):
 
diff --git a/sharpy/utils/solver_interface.py b/sharpy/utils/solver_interface.py
index eb671e496..898e4b70b 100644
--- a/sharpy/utils/solver_interface.py
+++ b/sharpy/utils/solver_interface.py
@@ -100,10 +100,8 @@ def dictionary_of_solvers(print_info=True):
     import sharpy.postproc
     dictionary = dict()
     for solver in dict_of_solvers:
-        if not solver.lower() == 'SaveData'.lower():
-            # TODO: why it does not work for savedata?
-            init_solver = initialise_solver(solver, print_info)
-            dictionary[solver] = init_solver.settings_default
+        init_solver = initialise_solver(solver, print_info)
+        dictionary[solver] = init_solver.settings_default
 
     return dictionary
 
@@ -193,5 +191,3 @@ def output_documentation(route=None):
                         out_file.write('    ./' + solver_type + '/' + k + '\n')
                 except AttributeError:
                     pass
-
-
diff --git a/utils/environment_linux.yml b/utils/environment_linux.yml
index f60c5bdc6..0834869ce 100644
--- a/utils/environment_linux.yml
+++ b/utils/environment_linux.yml
@@ -12,14 +12,15 @@ dependencies:
   - backcall=0.1.0=py37_0
   - bleach=3.1.0=py_0
   - bzip2=1.0.8=h7b6447c_0
-  - ca-certificates=2020.6.24=0
-  - certifi=2020.6.20=py37_0
+  - ca-certificates=2020.12.5=ha878542_0
+  - certifi=2020.12.5=py37h89c1867_1
   - cffi=1.13.1=py37h2e261b9_0
   - chardet=3.0.4=py37_1003
   - cmake=3.14.0=h52cb24c_0
   - colorama=0.4.1=py37_0
   - commonmark=0.9.0=py_0
   - configobj=5.0.6=py37_1
+  - control=0.8.4=py37h89c1867_0
   - cryptography=2.8=py37h1ba5d50_0
   - curl=7.67.0=hbc83047_0
   - cycler=0.10.0=py37_0
@@ -89,6 +90,7 @@ dependencies:
   - lz4-c=1.8.1.2=h14c3975_0
   - markupsafe=1.1.1=py37h7b6447c_0
   - matplotlib=3.1.1=py37h5429711_0
+  - matplotlib-base=3.1.3=py37hef1b27d_0
   - mayavi=4.7.1=py37h94891b3_2
   - mistune=0.8.4=py37h516909a_1000
   - more-itertools=7.2.0=py_0
@@ -98,7 +100,7 @@ dependencies:
   - ncurses=6.1=he6710b0_1
   - notebook=6.0.1=py37_0
   - numpy=1.17.3=py37h95a1406_0
-  - openssl=1.1.1g=h7b6447c_0
+  - openssl=1.1.1i=h27cfd23_0
   - packaging=19.2=py_0
   - pandas=0.25.3=py37hb3f55d8_0
   - pandoc=2.7.3=0
@@ -121,6 +123,7 @@ dependencies:
   - pysocks=1.7.1=py37_0
   - python=3.7.5=h0371630_0
   - python-dateutil=2.8.1=py_0
+  - python_abi=3.7=1_cp37m
   - pytz=2019.3=py_0
   - pyyaml=5.1.2=py37h7b6447c_0
   - pyzmq=18.1.0=py37he6710b0_0
@@ -136,13 +139,13 @@ dependencies:
   - six=1.13.0=py37_0
   - snowballstemmer=2.0.0=py_0
   - sphinx=3.0.3=py_0
+  - sphinx_rtd_theme=0.5.0=pyh9f0ad1d_0
   - sphinxcontrib-applehelp=1.0.1=py_0
   - sphinxcontrib-devhelp=1.0.1=py_0
   - sphinxcontrib-htmlhelp=1.0.2=py_0
   - sphinxcontrib-jsmath=1.0.1=py_0
   - sphinxcontrib-qthelp=1.0.2=py_0
   - sphinxcontrib-serializinghtml=1.1.3=py_0
-  - sphinx_rtd_theme=0.5.0=pyh9f0ad1d_0
   - sqlite=3.30.1=h7b6447c_0
   - tbb=2019.8=hfd86e86_0
   - terminado=0.8.3=py37_0
@@ -164,3 +167,4 @@ dependencies:
   - zipp=0.6.0=py_0
   - zlib=1.2.11=h7b6447c_3
   - zstd=1.3.7=h0b5b093_0
+  - slycot=0.4.0.0=py37h27181d0_1
diff --git a/utils/environment_minimal.yml b/utils/environment_minimal.yml
index 7d99ef79a..ea5841579 100644
--- a/utils/environment_minimal.yml
+++ b/utils/environment_minimal.yml
@@ -6,11 +6,12 @@ dependencies:
   - _libgcc_mutex=0.1=main
   - apptools=4.5.0=py_0
   - bzip2=1.0.8=h7b6447c_0
-  - ca-certificates=2020.6.24=0
-  - certifi=2020.6.20=py37_0
+  - ca-certificates=2020.12.5=ha878542_0
+  - certifi=2020.12.5=py37h89c1867_1
   - cmake=3.14.0=h52cb24c_0
   - colorama=0.4.1=py37_0
   - configobj=5.0.6=py_0
+  - control=0.8.4=py37h89c1867_0
   - curl=7.67.0=hbc83047_0
   - cycler=0.10.0=py_2
   - dbus=1.13.6=he372182_0
@@ -60,10 +61,11 @@ dependencies:
   - lxml=4.4.2=py37hefd8a0e_0
   - lz4-c=1.8.1.2=h14c3975_0
   - matplotlib=3.1.1=py37h5429711_0
+  - matplotlib-base=3.1.3=py37hef1b27d_0
   - mayavi=4.7.1=py37h94891b3_2
   - ncurses=6.1=he6710b0_1
   - numpy=1.17.3=py37h95a1406_0
-  - openssl=1.1.1g=h7b6447c_0
+  - openssl=1.1.1i=h27cfd23_0
   - pandas=0.25.3=py37hb3f55d8_0
   - pcre=8.43=he1b5a44_0
   - pip=19.3.1=py37_0
@@ -74,6 +76,7 @@ dependencies:
   - pyqt=5.9.2=py37hcca6a23_4
   - python=3.7.5=h0371630_0
   - python-dateutil=2.8.1=py_0
+  - python_abi=3.7=1_cp37m
   - pytz=2019.3=py_0
   - pyyaml=5.1.2=py37h7b6447c_0
   - qt=5.9.7=h52cfd70_2
@@ -83,6 +86,7 @@ dependencies:
   - setuptools=42.0.2=py37_0
   - sip=4.19.8=py37hf484d3e_0
   - six=1.13.0=py37_0
+  - slycot=0.4.0.0=py37h27181d0_1
   - sqlite=3.30.1=h7b6447c_0
   - tbb=2019.8=hfd86e86_0
   - tk=8.6.8=hbc83047_0
@@ -98,4 +102,3 @@ dependencies:
   - yaml=0.1.7=had09818_2
   - zlib=1.2.11=h7b6447c_3
   - zstd=1.3.7=h0b5b093_0
-


Test patch:
diff --git a/tests/coupled/multibody/floating_forces/MooringLineFD.dat b/tests/coupled/multibody/floating_forces/MooringLineFD.dat
new file mode 100644
index 000000000..d830fd110
--- /dev/null
+++ b/tests/coupled/multibody/floating_forces/MooringLineFD.dat
@@ -0,0 +1,509 @@
+This mooring line force-displacement relationship was generated by FAST (v6.10b-jmj, 19-Dec-2008)
+
+* LINE CHARACTERISTICS OF LINE NUMBER:   1 *         
+------------------------------------------------     
+  Line characteristics number:  1                    
+
+DISTANCE    TENSION    H.TENSION    SUSPL    TEN.ANCH
+  (m)         (kN)        (kN)       (m)       (kN)  
+
+ 653.00       174.6        0.1      250.0        0.1
+ 653.50       174.6        0.1      250.1        0.1
+ 654.00       174.7        0.2      250.2        0.2
+ 654.50       174.7        0.3      250.3        0.3
+ 655.00       174.8        0.3      250.4        0.3
+ 655.50       174.9        0.4      250.5        0.4
+ 656.00       174.9        0.5      250.6        0.5
+ 656.50       175.0        0.5      250.7        0.5
+ 657.00       175.1        0.6      250.8        0.6
+ 657.50       175.2        0.7      250.9        0.7
+ 658.00       175.3        0.8      251.1        0.8
+ 658.50       175.4        0.9      251.2        0.9
+ 659.00       175.4        1.0      251.3        1.0
+ 659.50       175.5        1.0      251.4        1.0
+ 660.00       175.6        1.1      251.6        1.1
+ 660.50       175.7        1.2      251.7        1.2
+ 661.00       175.8        1.3      251.8        1.3
+ 661.50       175.9        1.4      252.0        1.4
+ 662.00       176.0        1.5      252.1        1.5
+ 662.50       176.1        1.6      252.3        1.6
+ 663.00       176.2        1.7      252.4        1.7
+ 663.50       176.3        1.8      252.6        1.8
+ 664.00       176.4        1.9      252.7        1.9
+ 664.50       176.5        2.1      252.9        2.1
+ 665.00       176.7        2.2      253.0        2.2
+ 665.50       176.8        2.3      253.2        2.3
+ 666.00       176.9        2.4      253.4        2.4
+ 666.50       177.0        2.5      253.5        2.5
+ 667.00       177.1        2.6      253.7        2.6
+ 667.50       177.2        2.8      253.9        2.8
+ 668.00       177.4        2.9      254.0        2.9
+ 668.50       177.5        3.0      254.2        3.0
+ 669.00       177.6        3.1      254.4        3.1
+ 669.50       177.7        3.2      254.6        3.2
+ 670.00       177.9        3.4      254.7        3.4
+ 670.50       178.0        3.5      254.9        3.5
+ 671.00       178.1        3.6      255.1        3.6
+ 671.50       178.3        3.8      255.3        3.8
+ 672.00       178.4        3.9      255.5        3.9
+ 672.50       178.5        4.0      255.7        4.0
+ 673.00       178.7        4.2      255.9        4.2
+ 673.50       178.8        4.3      256.1        4.3
+ 674.00       178.9        4.5      256.3        4.5
+ 674.50       179.1        4.6      256.5        4.6
+ 675.00       179.2        4.8      256.7        4.8
+ 675.50       179.4        4.9      256.9        4.9
+ 676.00       179.5        5.1      257.1        5.1
+ 676.50       179.7        5.2      257.3        5.2
+ 677.00       179.8        5.4      257.5        5.4
+ 677.50       180.0        5.5      257.7        5.5
+ 678.00       180.2        5.7      257.9        5.7
+ 678.50       180.3        5.8      258.2        5.8
+ 679.00       180.5        6.0      258.4        6.0
+ 679.50       180.6        6.2      258.6        6.2
+ 680.00       180.8        6.3      258.8        6.3
+ 680.50       181.0        6.5      259.1        6.5
+ 681.00       181.1        6.6      259.3        6.6
+ 681.50       181.3        6.8      259.5        6.8
+ 682.00       181.5        7.0      259.8        7.0
+ 682.50       181.6        7.2      260.0        7.2
+ 683.00       181.8        7.3      260.2        7.3
+ 683.50       182.0        7.5      260.5        7.5
+ 684.00       182.2        7.7      260.7        7.7
+ 684.50       182.4        7.9      261.0        7.9
+ 685.00       182.5        8.1      261.2        8.1
+ 685.50       182.7        8.2      261.5        8.2
+ 686.00       182.9        8.4      261.7        8.4
+ 686.50       183.1        8.6      262.0        8.6
+ 687.00       183.3        8.8      262.2        8.8
+ 687.50       183.5        9.0      262.5        9.0
+ 688.00       183.7        9.2      262.8        9.2
+ 688.50       183.9        9.4      263.0        9.4
+ 689.00       184.1        9.6      263.3        9.6
+ 689.50       184.3        9.8      263.6        9.8
+ 690.00       184.5       10.0      263.9       10.0
+ 690.50       184.7       10.2      264.1       10.2
+ 691.00       184.9       10.4      264.4       10.4
+ 691.50       185.1       10.6      264.7       10.6
+ 692.00       185.3       10.8      265.0       10.8
+ 692.50       185.5       11.0      265.3       11.0
+ 693.00       185.7       11.2      265.6       11.2
+ 693.50       185.9       11.5      265.9       11.5
+ 694.00       186.2       11.7      266.1       11.7
+ 694.50       186.4       11.9      266.4       11.9
+ 695.00       186.6       12.1      266.7       12.1
+ 695.50       186.8       12.4      267.0       12.4
+ 696.00       187.1       12.6      267.4       12.6
+ 696.50       187.3       12.8      267.7       12.8
+ 697.00       187.5       13.0      268.0       13.0
+ 697.50       187.8       13.3      268.3       13.3
+ 698.00       188.0       13.5      268.6       13.5
+ 698.50       188.2       13.8      268.9       13.8
+ 699.00       188.5       14.0      269.2       14.0
+ 699.50       188.7       14.2      269.6       14.2
+ 700.00       189.0       14.5      269.9       14.5
+ 700.50       189.2       14.7      270.2       14.7
+ 701.00       189.5       15.0      270.6       15.0
+ 701.50       189.7       15.2      270.9       15.2
+ 702.00       190.0       15.5      271.2       15.5
+ 702.50       190.2       15.8      271.6       15.8
+ 703.00       190.5       16.0      271.9       16.0
+ 703.50       190.8       16.3      272.3       16.3
+ 704.00       191.0       16.5      272.6       16.5
+ 704.50       191.3       16.8      273.0       16.8
+ 705.00       191.6       17.1      273.3       17.1
+ 705.50       191.8       17.4      273.7       17.4
+ 706.00       192.1       17.6      274.0       17.6
+ 706.50       192.4       17.9      274.4       17.9
+ 707.00       192.7       18.2      274.8       18.2
+ 707.50       193.0       18.5      275.1       18.5
+ 708.00       193.2       18.8      275.5       18.8
+ 708.50       193.5       19.1      275.9       19.1
+ 709.00       193.8       19.4      276.3       19.4
+ 709.50       194.1       19.6      276.6       19.6
+ 710.00       194.4       19.9      277.0       19.9
+ 710.50       194.7       20.2      277.4       20.2
+ 711.00       195.0       20.5      277.8       20.5
+ 711.50       195.3       20.9      278.2       20.9
+ 712.00       195.6       21.2      278.6       21.2
+ 712.50       196.0       21.5      279.0       21.5
+ 713.00       196.3       21.8      279.4       21.8
+ 713.50       196.6       22.1      279.8       22.1
+ 714.00       196.9       22.4      280.2       22.4
+ 714.50       197.2       22.8      280.6       22.8
+ 715.00       197.6       23.1      281.1       23.1
+ 715.50       197.9       23.4      281.5       23.4
+ 716.00       198.2       23.8      281.9       23.8
+ 716.50       198.6       24.1      282.3       24.1
+ 717.00       198.9       24.4      282.8       24.4
+ 717.50       199.2       24.8      283.2       24.8
+ 718.00       199.6       25.1      283.6       25.1
+ 718.50       199.9       25.5      284.1       25.5
+ 719.00       200.3       25.8      284.5       25.8
+ 719.50       200.7       26.2      285.0       26.2
+ 720.00       201.0       26.5      285.4       26.5
+ 720.50       201.4       26.9      285.9       26.9
+ 721.00       201.7       27.3      286.3       27.3
+ 721.50       202.1       27.6      286.8       27.6
+ 722.00       202.5       28.0      287.3       28.0
+ 722.50       202.9       28.4      287.7       28.4
+ 723.00       203.2       28.8      288.2       28.8
+ 723.50       203.6       29.2      288.7       29.2
+ 724.00       204.0       29.5      289.2       29.5
+ 724.50       204.4       29.9      289.7       29.9
+ 725.00       204.8       30.3      290.1       30.3
+ 725.50       205.2       30.7      290.6       30.7
+ 726.00       205.6       31.1      291.1       31.1
+ 726.50       206.0       31.6      291.6       31.6
+ 727.00       206.4       32.0      292.1       32.0
+ 727.50       206.8       32.4      292.6       32.4
+ 728.00       207.3       32.8      293.2       32.8
+ 728.50       207.7       33.2      293.7       33.2
+ 729.00       208.1       33.6      294.2       33.6
+ 729.50       208.5       34.1      294.7       34.1
+ 730.00       209.0       34.5      295.3       34.5
+ 730.50       209.4       35.0      295.8       35.0
+ 731.00       209.9       35.4      296.3       35.4
+ 731.50       210.3       35.9      296.9       35.9
+ 732.00       210.8       36.3      297.4       36.3
+ 732.50       211.2       36.8      298.0       36.8
+ 733.00       211.7       37.2      298.5       37.2
+ 733.50       212.2       37.7      299.1       37.7
+ 734.00       212.6       38.2      299.6       38.2
+ 734.50       213.1       38.6      300.2       38.6
+ 735.00       213.6       39.1      300.8       39.1
+ 735.50       214.1       39.6      301.4       39.6
+ 736.00       214.6       40.1      301.9       40.1
+ 736.50       215.1       40.6      302.5       40.6
+ 737.00       215.6       41.1      303.1       41.1
+ 737.50       216.1       41.6      303.7       41.6
+ 738.00       216.6       42.1      304.3       42.1
+ 738.50       217.1       42.6      304.9       42.6
+ 739.00       217.6       43.1      305.5       43.1
+ 739.50       218.1       43.7      306.1       43.7
+ 740.00       218.7       44.2      306.8       44.2
+ 740.50       219.2       44.7      307.4       44.7
+ 741.00       219.7       45.3      308.0       45.3
+ 741.50       220.3       45.8      308.7       45.8
+ 742.00       220.8       46.4      309.3       46.4
+ 742.50       221.4       46.9      309.9       46.9
+ 743.00       222.0       47.5      310.6       47.5
+ 743.50       222.5       48.1      311.2       48.1
+ 744.00       223.1       48.6      311.9       48.6
+ 744.50       223.7       49.2      312.6       49.2
+ 745.00       224.3       49.8      313.2       49.8
+ 745.50       224.9       50.4      313.9       50.4
+ 746.00       225.5       51.0      314.6       51.0
+ 746.50       226.1       51.6      315.3       51.6
+ 747.00       226.7       52.2      316.0       52.2
+ 747.50       227.3       52.8      316.7       52.8
+ 748.00       227.9       53.5      317.4       53.5
+ 748.50       228.6       54.1      318.1       54.1
+ 749.00       229.2       54.7      318.8       54.7
+ 749.50       229.8       55.4      319.5       55.4
+ 750.00       230.5       56.0      320.3       56.0
+ 750.50       231.1       56.7      321.0       56.7
+ 751.00       231.8       57.3      321.7       57.3
+ 751.50       232.5       58.0      322.5       58.0
+ 752.00       233.2       58.7      323.2       58.7
+ 752.50       233.8       59.4      324.0       59.4
+ 753.00       234.5       60.1      324.8       60.1
+ 753.50       235.2       60.8      325.5       60.8
+ 754.00       235.9       61.5      326.3       61.5
+ 754.50       236.7       62.2      327.1       62.2
+ 755.00       237.4       62.9      327.9       62.9
+ 755.50       238.1       63.7      328.7       63.7
+ 756.00       238.8       64.4      329.5       64.4
+ 756.50       239.6       65.1      330.3       65.1
+ 757.00       240.3       65.9      331.1       65.9
+ 757.50       241.1       66.7      331.9       66.7
+ 758.00       241.9       67.4      332.7       67.4
+ 758.50       242.7       68.2      333.6       68.2
+ 759.00       243.4       69.0      334.4       69.0
+ 759.50       244.2       69.8      335.3       69.8
+ 760.00       245.0       70.6      336.1       70.6
+ 760.50       245.9       71.4      337.0       71.4
+ 761.00       246.7       72.2      337.9       72.2
+ 761.50       247.5       73.1      338.8       73.1
+ 762.00       248.4       73.9      339.6       73.9
+ 762.50       249.2       74.7      340.5       74.7
+ 763.00       250.1       75.6      341.4       75.6
+ 763.50       250.9       76.5      342.3       76.5
+ 764.00       251.8       77.4      343.3       77.4
+ 764.50       252.7       78.2      344.2       78.2
+ 765.00       253.6       79.1      345.1       79.1
+ 765.50       254.5       80.1      346.1       80.1
+ 766.00       255.4       81.0      347.0       81.0
+ 766.50       256.4       81.9      348.0       81.9
+ 767.00       257.3       82.8      348.9       82.8
+ 767.50       258.2       83.8      349.9       83.8
+ 768.00       259.2       84.8      350.9       84.8
+ 768.50       260.2       85.7      351.9       85.7
+ 769.00       261.2       86.7      352.9       86.7
+ 769.50       262.2       87.7      353.9       87.7
+ 770.00       263.2       88.7      354.9       88.7
+ 770.50       264.2       89.7      355.9       89.7
+ 771.00       265.2       90.8      357.0       90.8
+ 771.50       266.3       91.8      358.0       91.8
+ 772.00       267.3       92.9      359.1       92.9
+ 772.50       268.4       93.9      360.1       93.9
+ 773.00       269.5       95.0      361.2       95.0
+ 773.50       270.5       96.1      362.3       96.1
+ 774.00       271.7       97.2      363.4       97.2
+ 774.50       272.8       98.3      364.5       98.3
+ 775.00       273.9       99.5      365.6       99.5
+ 775.50       275.1      100.6      366.7      100.6
+ 776.00       276.2      101.8      367.8      101.8
+ 776.50       277.4      102.9      369.0      102.9
+ 777.00       278.6      104.1      370.1      104.1
+ 777.50       279.8      105.3      371.3      105.3
+ 778.00       281.0      106.6      372.4      106.6
+ 778.50       282.2      107.8      373.6      107.8
+ 779.00       283.5      109.0      374.8      109.0
+ 779.50       284.7      110.3      376.0      110.3
+ 780.00       286.0      111.6      377.2      111.6
+ 780.50       287.3      112.9      378.5      112.9
+ 781.00       288.6      114.2      379.7      114.2
+ 781.50       289.9      115.5      381.0      115.5
+ 782.00       291.3      116.9      382.2      116.9
+ 782.50       292.6      118.2      383.5      118.2
+ 783.00       294.0      119.6      384.8      119.6
+ 783.50       295.4      121.0      386.1      121.0
+ 784.00       296.8      122.4      387.4      122.4
+ 784.50       298.3      123.8      388.7      123.8
+ 785.00       299.7      125.3      390.0      125.3
+ 785.50       301.2      126.8      391.4      126.8
+ 786.00       302.7      128.2      392.7      128.2
+ 786.50       304.2      129.8      394.1      129.8
+ 787.00       305.7      131.3      395.5      131.3
+ 787.50       307.2      132.8      396.9      132.8
+ 788.00       308.8      134.4      398.3      134.4
+ 788.50       310.4      136.0      399.7      136.0
+ 789.00       312.0      137.6      401.1      137.6
+ 789.50       313.6      139.2      402.6      139.2
+ 790.00       315.3      140.9      404.1      140.9
+ 790.50       317.0      142.5      405.5      142.5
+ 791.00       318.7      144.2      407.0      144.2
+ 791.50       320.4      146.0      408.5      146.0
+ 792.00       322.1      147.7      410.1      147.7
+ 792.50       323.9      149.5      411.6      149.5
+ 793.00       325.7      151.3      413.2      151.3
+ 793.50       327.5      153.1      414.7      153.1
+ 794.00       329.3      154.9      416.3      154.9
+ 794.50       331.2      156.8      417.9      156.8
+ 795.00       333.1      158.7      419.5      158.7
+ 795.50       335.0      160.6      421.2      160.6
+ 796.00       336.9      162.5      422.8      162.5
+ 796.50       338.9      164.5      424.5      164.5
+ 797.00       340.9      166.5      426.1      166.5
+ 797.50       342.9      168.5      427.8      168.5
+ 798.00       345.0      170.6      429.6      170.6
+ 798.50       347.1      172.7      431.3      172.7
+ 799.00       349.2      174.8      433.0      174.8
+ 799.50       351.4      176.9      434.8      176.9
+ 800.00       353.5      179.1      436.6      179.1
+ 800.50       355.7      181.3      438.4      181.3
+ 801.00       358.0      183.6      440.2      183.6
+ 801.50       360.3      185.9      442.1      185.9
+ 802.00       362.6      188.2      443.9      188.2
+ 802.50       364.9      190.5      445.8      190.5
+ 803.00       367.3      192.9      447.7      192.9
+ 803.50       369.7      195.3      449.6      195.3
+ 804.00       372.1      197.7      451.6      197.7
+ 804.50       374.6      200.2      453.6      200.2
+ 805.00       377.1      202.8      455.5      202.8
+ 805.50       379.7      205.3      457.5      205.3
+ 806.00       382.3      207.9      459.6      207.9
+ 806.50       385.0      210.6      461.6      210.6
+ 807.00       387.6      213.2      463.7      213.2
+ 807.50       390.4      216.0      465.8      216.0
+ 808.00       393.1      218.7      467.9      218.7
+ 808.50       395.9      221.5      470.1      221.5
+ 809.00       398.8      224.4      472.2      224.4
+ 809.50       401.7      227.3      474.4      227.3
+ 810.00       404.6      230.2      476.6      230.2
+ 810.50       407.6      233.2      478.9      233.2
+ 811.00       410.6      236.3      481.1      236.3
+ 811.50       413.7      239.4      483.4      239.4
+ 812.00       416.9      242.5      485.7      242.5
+ 812.50       420.1      245.7      488.1      245.7
+ 813.00       423.3      248.9      490.4      248.9
+ 813.50       426.6      252.2      492.8      252.2
+ 814.00       429.9      255.6      495.2      255.6
+ 814.50       433.3      259.0      497.7      259.0
+ 815.00       436.8      262.4      500.2      262.4
+ 815.50       440.3      265.9      502.7      265.9
+ 816.00       443.9      269.5      505.2      269.5
+ 816.50       447.5      273.2      507.8      273.2
+ 817.00       451.2      276.9      510.4      276.9
+ 817.50       455.0      280.6      513.0      280.6
+ 818.00       458.8      284.4      515.7      284.4
+ 818.50       462.7      288.3      518.4      288.3
+ 819.00       466.6      292.3      521.1      292.3
+ 819.50       470.7      296.3      523.8      296.3
+ 820.00       474.8      300.4      526.6      300.4
+ 820.50       478.9      304.6      529.4      304.6
+ 821.00       483.2      308.8      532.3      308.8
+ 821.50       487.5      313.2      535.2      313.2
+ 822.00       491.9      317.6      538.1      317.6
+ 822.50       496.4      322.0      541.1      322.0
+ 823.00       500.9      326.6      544.1      326.6
+ 823.50       505.6      331.2      547.1      331.2
+ 824.00       510.3      335.9      550.2      335.9
+ 824.50       515.1      340.8      553.3      340.8
+ 825.00       520.0      345.7      556.5      345.7
+ 825.50       525.0      350.6      559.7      350.6
+ 826.00       530.0      355.7      562.9      355.7
+ 826.50       535.2      360.9      566.2      360.9
+ 827.00       540.5      366.2      569.5      366.2
+ 827.50       545.9      371.5      572.8      371.5
+ 828.00       551.3      377.0      576.2      377.0
+ 828.50       556.9      382.6      579.7      382.6
+ 829.00       562.6      388.3      583.2      388.3
+ 829.50       568.4      394.1      586.7      394.1
+ 830.00       574.3      400.0      590.3      400.0
+ 830.50       580.3      406.0      593.9      406.0
+ 831.00       586.5      412.2      597.6      412.2
+ 831.50       592.7      418.4      601.4      418.4
+ 832.00       599.1      424.8      605.1      424.8
+ 832.50       605.6      431.3      609.0      431.3
+ 833.00       612.3      438.0      612.9      438.0
+ 833.50       619.0      444.8      616.8      444.8
+ 834.00       626.0      451.7      620.8      451.7
+ 834.50       633.0      458.7      624.8      458.7
+ 835.00       640.2      465.9      628.9      465.9
+ 835.50       647.6      473.3      633.1      473.3
+ 836.00       655.1      480.8      637.3      480.8
+ 836.50       662.7      488.5      641.6      488.5
+ 837.00       670.6      496.3      645.9      496.3
+ 837.50       678.5      504.3      650.3      504.3
+ 838.00       686.7      512.4      654.8      512.4
+ 838.50       695.0      520.8      659.3      520.8
+ 839.00       703.5      529.3      663.9      529.3
+ 839.50       712.2      538.0      668.6      538.0
+ 840.00       721.1      546.9      673.3      546.9
+ 840.50       730.2      555.9      678.1      555.9
+ 841.00       739.4      565.2      683.0      565.2
+ 841.50       748.9      574.7      687.9      574.7
+ 842.00       758.6      584.4      692.9      584.4
+ 842.50       768.5      594.3      698.0      594.3
+ 843.00       778.6      604.4      703.1      604.4
+ 843.50       789.0      614.8      708.4      614.8
+ 844.00       799.6      625.4      713.7      625.4
+ 844.50       810.4      636.2      719.1      636.2
+ 845.00       821.5      647.3      724.6      647.3
+ 845.50       832.8      658.6      730.1      658.6
+ 846.00       844.4      670.2      735.8      670.2
+ 846.50       856.3      682.1      741.5      682.1
+ 847.00       868.4      694.2      747.3      694.2
+ 847.50       880.8      706.7      753.2      706.7
+ 848.00       893.6      719.4      759.2      719.4
+ 848.50       906.6      732.4      765.3      732.4
+ 849.00       919.9      745.8      771.5      745.8
+ 849.50       933.6      759.4      777.8      759.4
+ 850.00       947.6      773.4      784.2      773.4
+ 850.50       961.9      787.8      790.7      787.8
+ 851.00       976.6      802.5      797.3      802.5
+ 851.50       991.6      817.5      804.0      817.5
+ 852.00      1007.0      832.9      810.8      832.9
+ 852.50      1022.8      848.7      817.7      848.7
+ 853.00      1039.0      864.9      824.7      864.9
+ 853.50      1055.6      881.5      831.8      881.5
+ 854.00      1072.6      898.5      839.1      898.5
+ 854.50      1090.0      916.0      846.5      916.0
+ 855.00      1107.9      933.9      854.0      933.9
+ 855.50      1126.3      952.2      861.6      952.2
+ 856.00      1145.1      971.0      869.3      971.0
+ 856.50      1164.4      990.3      877.2      990.3
+ 857.00      1184.1     1010.1      885.2     1010.1
+ 857.50      1204.4     1030.4      893.3     1030.4
+ 858.00      1225.3     1051.3      901.6     1051.3
+ 858.50      1246.9     1072.9      902.2     1072.9
+ 859.00      1269.7     1095.7      902.2     1095.8
+ 859.50      1293.9     1119.7      902.2     1119.9
+ 860.00      1319.4     1145.2      902.2     1145.5
+ 860.50      1346.5     1172.1      902.2     1172.6
+ 861.00      1375.3     1200.7      902.2     1201.4
+ 861.50      1405.9     1231.0      902.2     1232.0
+ 862.00      1438.5     1263.2      902.2     1264.6
+ 862.50      1473.3     1297.6      902.2     1299.4
+ 863.00      1510.5     1334.3      902.2     1336.6
+ 863.50      1550.3     1373.5      902.2     1376.4
+ 864.00      1593.0     1415.6      902.2     1419.1
+ 864.50      1638.8     1460.6      902.2     1465.0
+ 865.00      1688.1     1509.0      902.2     1514.3
+ 865.50      1741.2     1561.0      902.2     1567.4
+ 866.00      1798.4     1617.0      902.2     1624.6
+ 866.50      1860.0     1677.3      902.2     1686.3
+ 867.00      1926.6     1742.3      902.2     1752.9
+ 867.50      1998.4     1812.4      902.2     1824.8
+ 868.00      2075.9     1887.9      902.2     1902.3
+ 868.50      2159.4     1969.2      902.2     1985.8
+ 869.00      2249.2     2056.5      902.2     2075.6
+ 869.50      2345.7     2150.3      902.2     2172.2
+ 870.00      2449.0     2250.6      902.2     2275.5
+ 870.50      2559.3     2357.7      902.2     2385.9
+ 871.00      2676.6     2471.5      902.2     2503.3
+ 871.50      2801.0     2591.9      902.2     2627.7
+ 872.00      2932.1     2718.9      902.2     2758.9
+ 872.50      3069.8     2852.2      902.2     2896.7
+ 873.00      3213.8     2991.5      902.2     3040.7
+ 873.50      3363.6     3136.4      902.2     3190.6
+ 874.00      3518.9     3286.5      902.2     3345.9
+ 874.50      3679.2     3441.4      902.2     3506.3
+ 875.00      3844.0     3600.6      902.2     3671.2
+ 875.50      4012.9     3763.8      902.2     3840.2
+ 876.00      4185.6     3930.5      902.2     4012.9
+ 876.50      4361.5     4100.4      902.2     4188.9
+ 877.00      4540.4     4273.0      902.2     4367.9
+ 877.50      4721.9     4448.2      902.2     4549.5
+ 878.00      4905.7     4625.7      902.2     4733.4
+ 878.50      5091.7     4805.1      902.2     4919.4
+ 879.00      5279.4     4986.2      902.2     5107.2
+ 879.50      5468.8     5169.0      902.2     5296.7
+ 880.00      5659.6     5353.1      902.2     5487.6
+ 880.50      5851.7     5538.4      902.2     5679.7
+ 881.00      6044.9     5724.9      902.2     5873.1
+ 881.50      6239.2     5912.3      902.2     6067.4
+ 882.00      6434.4     6100.6      902.2     6262.7
+ 882.50      6630.4     6289.7      902.2     6458.8
+ 883.00      6827.1     6479.5      902.2     6655.6
+ 883.50      7024.5     6670.0      902.2     6853.1
+ 884.00      7222.5     6861.0      902.2     7051.2
+ 884.50      7421.0     7052.5      902.2     7249.8
+ 885.00      7620.1     7244.6      902.2     7448.9
+ 885.50      7819.5     7437.0      902.2     7648.4
+ 886.00      8019.4     7629.9      902.2     7848.4
+ 886.50      8219.6     7823.1      902.2     8048.7
+ 887.00      8420.1     8016.6      902.2     8249.3
+ 887.50      8621.0     8210.4      902.2     8450.3
+ 888.00      8822.1     8404.5      902.2     8651.5
+ 888.50      9023.5     8598.9      902.2     8853.0
+ 889.00      9225.1     8793.5      902.2     9054.7
+ 889.50      9427.0     8988.3      902.2     9256.6
+ 890.00      9629.0     9183.4      902.2     9458.7
+ 890.50      9831.2     9378.6      902.2     9661.0
+ 891.00     10033.6     9574.0      902.2     9863.5
+ 891.50     10236.2     9769.5      902.2    10066.1
+ 892.00     10438.9     9965.2      902.2    10268.9
+ 892.50     10641.7    10161.1      902.2    10471.8
+ 893.00     10844.7    10357.1      902.2    10674.9
+ 893.50     11047.7    10553.2      902.2    10878.1
+ 894.00     11250.9    10749.4      902.2    11081.3
+ 894.50     11454.2    10945.8      902.2    11284.7
+ 895.00     11657.6    11142.2      902.2    11488.2
+ 895.50     11861.1    11338.8      902.2    11691.8
+ 896.00     12064.7    11535.4      902.2    11895.4
+ 896.50     12268.3    11732.1      902.2    12099.2
+ 897.00     12472.0    11929.0      902.2    12303.0
+ 897.50     12675.8    12125.9      902.2    12506.9
+ 898.00     12879.7    12322.9      902.2    12710.8
+ 898.50     13083.6    12519.9      902.2    12914.8
+ 899.00     13287.6    12717.1      902.2    13118.9
+ 899.50     13491.7    12914.3      902.2    13323.1
+ 900.00     13695.8    13111.5      902.2    13527.2
+ 900.50     13900.0    13308.9      902.2    13731.5
+ 901.00     14104.2    13506.2      902.2    13935.8
+ 901.50     14308.4    13703.7      902.2    14140.1
+ 902.00     14512.7    13901.2      902.2    14344.5
+ 902.50     14717.1    14098.8      902.2    14548.9
diff --git a/tests/coupled/multibody/floating_forces/test_floatingforces.py b/tests/coupled/multibody/floating_forces/test_floatingforces.py
new file mode 100644
index 000000000..fb8470c05
--- /dev/null
+++ b/tests/coupled/multibody/floating_forces/test_floatingforces.py
@@ -0,0 +1,176 @@
+import numpy as np
+import unittest
+import os
+import shutil
+from scipy import fft
+import sharpy.generators.floatingforces as ff
+
+
+class TestFloatingForces(unittest.TestCase):
+    """
+    Some tests of the floating forces library
+    Check references therein
+    """
+
+    def test_compute_xf_zf(self):
+        """
+            This function tests based on by hand computations and data from the MooringLineFD.txt file
+            from the OC3 task.
+                Jonkman, J.
+                Definition of the Floating System for Phase IV of OC3
+                NREL/TP-500-47535
+        """
+    
+        # Values for OC3
+        l = 902.2 # Initial length [m]
+        w = 698.094 # Aparent weight 77.7066*9.81 # Apparent mass per unit length times gravity
+        EA = 384243000. # Extensional stiffness
+        cb = 0.1 # Seabed friction coefficient
+    
+        # No mooring line on the Seabed
+        vf = 1.1*l*w # 692802.4475
+        hf = vf
+        xf_byhand = 784.5965853 + 1.626695524
+        zf_byhand = 406.9813526 + 0.887288467
+        xf, zf = ff.compute_xf_zf(hf, vf, l, w, EA, cb)
+        self.assertAlmostEqual(xf_byhand, xf, 4)
+        self.assertAlmostEqual(zf_byhand, zf, 4)
+    
+        # Some mooring line on the Seabed
+        lb_div_l = 0.1 # 10% of the mooring line on the seabed
+        vf = (1-lb_div_l)*l*w
+        hf = vf
+        xf_byhand = 90.22 + 715.6577252 + 1.330932701 - 7.298744381e-4
+        zf_byhand = 336.3331284 + 0.598919715
+        xf, zf = ff.compute_xf_zf(hf, vf, l, w, EA, cb)
+        self.assertAlmostEqual(xf_byhand, xf, 4)
+        self.assertAlmostEqual(zf_byhand, zf, 4)
+    
+    
+    def test_generate_mooringlinefd(self):
+        """
+            This function generates a file similar to MoorinLinesFD.txt for compiarison
+            File obtained from the google drive directory of the OC3 benchmark
+        """
+    
+        # Values for OC3
+        l = 902.2 # Initial length [m]
+        w = 698.094 # Aparent weight 77.7066*9.81 # Apparent mass per unit length times gravity
+        EA = 384243000. # Extensional stiffness
+        cb = 0. # Seabed friction coefficient
+    
+        zf = 320. - 70.
+    
+        # xf0 = 853.87
+        # xf_list = np.arange(653.00, 902.50 + 1., 1.)
+        xf_list = np.arange(653.00, 902.50, 10.)
+        npoints = xf_list.shape[0]
+        output = np.zeros((npoints, 4))
+        for i in range(npoints):
+            hf, vf = ff.quasisteady_mooring(xf_list[i], zf, l, w, EA, cb, hf0=None, vf0=None)
+            # print(xf0, zf0, hf0, vf0)
+            lb = np.maximum(l - vf/w, 0)
+            # print("Suspended lenght = %f" % (l - lb))
+            output[i, :] = np.array([xf_list[i], np.sqrt(vf**2 + hf**2)*1e-3, hf*1e-3, (l - lb)])
+
+        of_results = np.loadtxt("MooringLineFD.dat", skiprows=9)
+        for i in range(npoints):
+            for icol in range(4):
+                of_value = np.interp(xf_list[i], of_results[:, 0], of_results[:, icol])
+                error = np.abs((np.round_(output[i, icol], 1) - of_value)/of_value)
+                self.assertLess(error, 0.1)
+                
+        save_txt = False
+        if save_txt:    
+            np.savetxt("sharpy_mooringlinefd.txt", output, header="# DISTANCE(m) TENSION(kN) HTENSION(kN) SUSPL(m)")
+
+    
+    def test_change_system(self):
+        # Wind turbine degrees of freedom: Surge, sway, heave, roll, pitch, yaw.
+        # SHARPy axis associated:              z,    y,     x,    z,     y,   x
+    
+        wt_matrix_num = np.zeros((6,6),)
+        for idof in range(6):
+            for jdof in range(6):
+                wt_matrix_num[idof, jdof] = 10.*idof + jdof
+    
+        sharpy_matrix = ff.change_of_to_sharpy(wt_matrix_num)
+        undo_sharpy_matrix = ff.change_of_to_sharpy(sharpy_matrix)
+        
+        for idof in range(6):
+            for jdof in range(6):
+                self.assertEqual(wt_matrix_num[idof, jdof], undo_sharpy_matrix[idof, jdof])
+   
+ 
+    def test_time_wave_forces(self):
+        Tp = 14.656 #10.
+        Hs = 5.49 #6.
+        nrealisations = 100
+        dt = 2*np.pi/4 # 1./20 To get max freq equal to 10Hz
+        ntime_steps = 1000
+        time = np.arange(ntime_steps)*dt
+    
+        # Get the zero-noise specturm
+        w_js = np.arange(0, 4, 0.01)
+        zero_noise_spectrum = ff.jonswap_spectrum(Tp, Hs, w_js)
+    
+        # Compute different realisations
+        xi = np.zeros((2, 6), dtype=complex)
+        xi[0, 0] = 1. + 0j
+        xi[1, 0] = 1. + 0j
+        w_xi = np.array([0., 4.])
+        wave_force = np.zeros((ntime_steps, nrealisations), dtype=np.complex)
+        for ireal in range(nrealisations):
+            wave_force[:, ireal] = ff.time_wave_forces(Tp, Hs, dt, time, xi, w_xi)[:, 0] # Keep only on dimension
+    
+        # Compute the spectrum of the realisations
+        ns = np.zeros((ntime_steps//2, nrealisations), dtype=np.complex)
+        for ireal in range(nrealisations):
+            ns[:, ireal] = dt/ntime_steps*np.abs(fft(wave_force[:, ireal])[:ntime_steps//2])**2
+            ns[1:, ireal] *= 2
+        # To rad/s
+        ns /= 2.*np.pi
+        w_ns = (np.fft.fftfreq(ntime_steps, d=dt)[:ntime_steps//2])*2.*np.pi
+        # Compare the zero noise with the realisations average
+        avg_noise_spectrum = np.average(ns, axis=1)
+
+        for iomega in range(w_js.shape[0]):
+            error = (np.interp(w_js[iomega], w_ns, avg_noise_spectrum) - zero_noise_spectrum[iomega])
+            if error > 0.1:
+                error /= zero_noise_spectrum[iomega]
+                # 0.3 is a large error but otherwise I need to increment nrealisations a lot
+                # Use ``save_fig = True`` for visual inspection
+                self.assertLess(error, 0.3)
+            
+        save_fig = False
+        if save_fig:  
+            import matplotlib.pyplot as plt
+            fig, ax = plt.subplots(1, 1, figsize=(4, 3))
+            ax.grid()
+            ax.set_xlabel("omega [rad/s]")
+            ax.set_ylabel("specturm")
+            ax.set_xlim(0, 4)
+            ax.set_ylim(0, 12)
+            for ireal in range(nrealisations):
+                ax.plot(w_ns, np.abs(ns[:, ireal]), 'bo')
+            ax.plot(w_ns, avg_noise_spectrum, '-', label="avg")
+            ax.plot(w_js, zero_noise_spectrum, '--', label="JONSWAP")
+            fig.legend()
+            fig.tight_layout()
+            fig.show()
+            fig.savefig("spectrum.png")
+            plt.close()
+
+
+    # def tearDown(self):
+    #     solver_path = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))
+    #     solver_path += '/'
+    #     files_to_delete = [name + '.aero.h5',
+    #                        name + '.dyn.h5',
+    #                        name + '.fem.h5',
+    #                        name + '.mb.h5',
+    #                        name + '.sharpy']
+    #     for f in files_to_delete:
+    #         os.remove(solver_path + f)
+
+    #     shutil.rmtree(solver_path + 'output/')