diff --git a/README.md b/README.md
index 5f7ab4a19..75c434a80 100644
--- a/README.md
+++ b/README.md
@@ -33,6 +33,8 @@ vortex ring lattice with the boundary conditions enforced at the collocation poi
 The Kutta condition is also enforced at the trailing edge. The wake can be simulated by either additional vortex rings
 or by infinitely long horseshoe vortices, which are ideally suited for steady simulations only.
 
+The aerodynamic model has recently been extended by a linear source panel method (SPM) to model nonlifting bodies for example fuselages. The SPM and UVLM can be coupled to model fuselage-wing configuration and a junction handling approach, based on phantom panels and circulation interpolation, has been added.
+
 The input problems can be structural, aerodynamic or coupled, yielding an aeroelastic system.
 
 ## [Capabilities](http://ic-sharpy.readthedocs.io/en/latest/content/capabilities.html)
@@ -42,7 +44,7 @@ wings and wind turbines. In addition, it supports linearisation of these nonline
 arbitrary conditions and includes various tools such as: model reduction or frequency analysis.
 
 In short, SHARPy offers (amongst others) the following solutions to the user:
-* Static aerodynamic, structural and aeroelastic solutions
+* Static aerodynamic, structural and aeroelastic solutions including fuselage effects
 * Finding trim conditions for aeroelastic configurations
 * Nonlinear, dynamic time domain simulations under a large number of conditions such as:
     + Prescribed trajectories.
diff --git a/docs/source/content/casefiles.rst b/docs/source/content/casefiles.rst
index db537ac60..da570fc42 100644
--- a/docs/source/content/casefiles.rst
+++ b/docs/source/content/casefiles.rst
@@ -321,6 +321,55 @@ Item by item:
     should be included for each airfoil defined. Each entry consists of a 4-column table. The first column corresponds
     to the angle of attack (in radians) and then the ``C_L``, ``C_D`` and ``C_M``. 
 
+Nonlifting Body file
+-----------------
+
+All the nonlifting body data is contained in ``case.nonlifting_body.h5``.
+
+The idea behind the structure of the model definition of nonlifting bodies in SHARPy is similiar to the aerodynamic 
+one for lifting surfaces. Again for each node or element we define several parameters.
+
+Item by item:
+
+* ``shape``: Type of geometrical form of 3D nonlifting body.
+
+    In the ``nonlifting_body.h5`` file, there is a Group called ``shape``. The shape indicates the geometrical form of the 
+    nonlifting body. Common options for this parameter are ``'cylindrical'`` and ``'specific'``. For the former, SHARPy 
+    expects rotational symmetric cross-section for which only a radius is required for each node. For the ``'specific'`` 
+    option, SHARPy can create a more unique nonlifting body geometry by creating an ellipse at each fuselage defined by 
+    :math:`\frac{y^2}{a^2}+\frac{z^2}{b^2}=1` with the given ellipse axis lengths :math:`a`  and :math:`b`. Further, SHARPy 
+    lets define the user to create a vertical offset from the node with :math:`z_0`.
+
+* ``radius [num_node]``: Cross-sectional radius.
+
+    Is an array with the radius of specified for each fuselage node.
+
+* ``a_ellipse [num_node]``: Elliptical axis lengths along the local y-axis.
+
+    Is an array with the length of the elliptical axis along the y-axis.
+
+* ``b_ellipse [num_node]``: Elliptical axis lengths along the local z-axis.
+
+    Is an array with the length of the elliptical axis along the z-axis.
+
+* ``z_0_ellipse [num_node]``: Vertical offset of the ellipse center from the beam node.
+
+    Is an array with the vertical offset of the center of the elliptical cross-sectoin from the fuselage node.
+
+*  ``surface_m [num_surfaces]``: Radial panelling.
+
+    Is an integer array with the number of radial panels for every surface.
+
+*  ``nonlifting_body_node [num_node]``: Nonlifting body node definition.
+
+    Is a boolean (``True`` or ``False``) array that indicates if that node has a nonlifting body
+    attached to it.
+    
+*  ``surface_distribution [num_elem]``:  Nonlifting Surface integer array.
+
+    It contains the index of the surface the element belongs to. Surfaces need to be continuous, so please note 
+    that if your beam numbering is not continuous, you need to make a surface per continuous section.
+
 
 Time-varying force input file (``.dyn.h5``)
 -------------------------------------------
diff --git a/docs/source/content/example_notebooks/wind_turbine.ipynb b/docs/source/content/example_notebooks/wind_turbine.ipynb
index be366343a..3d1722c5a 100644
--- a/docs/source/content/example_notebooks/wind_turbine.ipynb
+++ b/docs/source/content/example_notebooks/wind_turbine.ipynb
@@ -1010,7 +1010,7 @@
     "        inode_in_elem = sharpy_output.structure.node_master_elem[node_global_index, 1]\n",
     "        CAB = algebra.crv2rotation(tstep.psi[ielem, inode_in_elem, :])\n",
     "\n",
-    "        c[iblade][inode] = sharpy_output.aero.aero_dict['chord'][ielem,inode_in_elem]\n",
+    "        c[iblade][inode] = sharpy_output.aero.data_dict['chord'][ielem,inode_in_elem]\n",
     "\n",
     "        forces_AFoR = np.dot(CAB, forces[iblade][inode, 0:3])\n",
     "\n",
diff --git a/docs/source/index.rst b/docs/source/index.rst
index 6ae696e07..adbc7cf5d 100644
--- a/docs/source/index.rst
+++ b/docs/source/index.rst
@@ -29,7 +29,7 @@ SHARPy is an aeroelastic analysis package currently under development at the Dep
 Imperial College London. It can be used for the structural, aerodynamic, aeroelastic and flight dynamics
 analysis of flexible aircraft, flying wings and wind turbines. Amongst other capabilities_, it offers the following solutions to the user:
 
-* Static aerodynamic, structural and aeroelastic solutions
+* Static aerodynamic, structural and aeroelastic solutions including fuselage effects
 * Finding trim conditions for aeroelastic configurations
 * Nonlinear, dynamic time domain simulations under a large number of conditions such as:
 
diff --git a/lib/UVLM b/lib/UVLM
index 23c13b6d4..d8af34a22 160000
--- a/lib/UVLM
+++ b/lib/UVLM
@@ -1 +1 @@
-Subproject commit 23c13b6d4d7da247b97b34ee06d29172c159fb42
+Subproject commit d8af34a22baf1cddd38f1e362274c407637aab1c
diff --git a/sharpy/aero/models/aerogrid.py b/sharpy/aero/models/aerogrid.py
index 399ce6c0c..4ba0d4f74 100644
--- a/sharpy/aero/models/aerogrid.py
+++ b/sharpy/aero/models/aerogrid.py
@@ -14,8 +14,10 @@
 from sharpy.utils.datastructures import AeroTimeStepInfo
 import sharpy.utils.generator_interface as gen_interface
 
+from sharpy.aero.models.grid import Grid
 
-class Aerogrid(object):
+
+class Aerogrid(Grid):
     """
     ``Aerogrid`` is the main object containing information of the grid of panels
 
@@ -23,70 +25,34 @@ class Aerogrid(object):
 
     """
     def __init__(self):
-        self.aero_dict = None
-        self.beam = None
-        self.aero_settings = None
-
-        self.timestep_info = []
-        self.ini_info = None
-
-        self.surface_distribution = None
-        self.surface_m = None
-        self.aero_dimensions = None
-        self.aero_dimensions_star = None
+        super().__init__()
+        self.dimensions_star = None
         self.airfoil_db = dict()
-        self.struct2aero_mapping = None
-        self.aero2struct_mapping = []
 
-        self.n_node = 0
-        self.n_elem = 0
-        self.n_surf = 0
-        self.n_aero_node = 0
+        self.grid_type = "aero"
         self.n_control_surfaces = 0
 
         self.cs_generators = []
 
-        self.polars = None
-        self.wake_shape_generator = None
-
-    def generate(self, aero_dict, beam, aero_settings, ts):
-        self.aero_dict = aero_dict
-        self.beam = beam
-        self.aero_settings = aero_settings
-
-        # number of total nodes (structural + aero&struc)
-        self.n_node = len(aero_dict['aero_node'])
-        # number of elements
-        self.n_elem = len(aero_dict['surface_distribution'])
-        # surface distribution
-        self.surface_distribution = aero_dict['surface_distribution']
-        # number of surfaces
-        temp = set(aero_dict['surface_distribution'])
-        self.n_surf = sum(1 for i in temp if i >= 0)
-        # number of chordwise panels
-        self.surface_m = aero_dict['surface_m']
-        # number of aero nodes
-        self.n_aero_node = sum(aero_dict['aero_node'])
-
-        # get N per surface
-        self.calculate_dimensions()
+    def generate(self, data_dict, beam, settings, ts):
+        super().generate(data_dict, beam, settings, ts)
 
         # write grid info to screen
         self.output_info()
 
         # allocating initial grid storage
-        self.ini_info = AeroTimeStepInfo(self.aero_dimensions,
-                                         self.aero_dimensions_star)
+        self.ini_info = AeroTimeStepInfo(self.dimensions,
+                                         self.dimensions_star)
 
         # load airfoils db
         # for i_node in range(self.n_node):
         for i_elem in range(self.n_elem):
             for i_local_node in range(self.beam.num_node_elem):
                 try:
-                    self.airfoil_db[self.aero_dict['airfoil_distribution'][i_elem, i_local_node]]
+                    self.airfoil_db[self.data_dict['airfoil_distribution'][i_elem, i_local_node]]
                 except KeyError:
-                    airfoil_coords = self.aero_dict['airfoils'][str(self.aero_dict['airfoil_distribution'][i_elem, i_local_node])]
-                    self.airfoil_db[self.aero_dict['airfoil_distribution'][i_elem, i_local_node]] = (
+                    airfoil_coords = self.data_dict['airfoils'][str(self.data_dict['airfoil_distribution'][i_elem, i_local_node])]
+                    self.airfoil_db[self.data_dict['airfoil_distribution'][i_elem, i_local_node]] = (
                         scipy.interpolate.interp1d(airfoil_coords[:, 0],
                                                    airfoil_coords[:, 1],
                                                    kind='quadratic',
@@ -94,39 +60,40 @@ def generate(self, aero_dict, beam, aero_settings, ts):
                                                    fill_value='extrapolate',
                                                    assume_sorted=True))
         try:
-            self.n_control_surfaces = np.sum(np.unique(self.aero_dict['control_surface']) >= 0)
+            self.n_control_surfaces = np.sum(np.unique(self.data_dict['control_surface']) >= 0)
         except KeyError:
             pass
 
-        # Backward compatibility: check whether control surface deflection settings have been specified. If not, create
+        # Backward compatibility: check whether control surface deflection aero_settings have been specified. If not, create
         # section with empty list, such that no cs generator is appended
         try:
-            aero_settings['control_surface_deflection']
+            settings['control_surface_deflection']
         except KeyError:
-            aero_settings.update({'control_surface_deflection': ['']*self.n_control_surfaces})
+            settings.update({'control_surface_deflection': ['']*self.n_control_surfaces})
 
         # pad ctrl surfaces dict with empty strings if not defined
-        if len(aero_settings['control_surface_deflection']) != self.n_control_surfaces:
-            undef_ctrl_sfcs = ['']*(self.n_control_surfaces - len(aero_settings['control_surface_deflection']))
-            aero_settings['control_surface_deflection'].extend(undef_ctrl_sfcs)
+        if len(settings['control_surface_deflection']) != self.n_control_surfaces:
+            undef_ctrl_sfcs = ['']*(self.n_control_surfaces - len(settings['control_surface_deflection']))
+            settings['control_surface_deflection'].extend(undef_ctrl_sfcs)
+
 
         # initialise generators
         with_error_initialising_cs = False
         for i_cs in range(self.n_control_surfaces):
-            if aero_settings['control_surface_deflection'][i_cs] == '':
+            if settings['control_surface_deflection'][i_cs] == '':
                 self.cs_generators.append(None)
             else:
                 cout.cout_wrap('Initialising Control Surface {:g} generator'.format(i_cs), 1)
                 # check that the control surface is not static
-                if self.aero_dict['control_surface_type'][i_cs] == 0:
+                if self.data_dict['control_surface_type'][i_cs] == 0:
                     raise TypeError('Control surface {:g} is defined as static but there is a control surface generator'
                                     'associated with it'.format(i_cs))
                 generator_type = gen_interface.generator_from_string(
-                    aero_settings['control_surface_deflection'][i_cs])
+                    settings['control_surface_deflection'][i_cs])
                 self.cs_generators.append(generator_type())
                 try:
                     self.cs_generators[i_cs].initialise(
-                        aero_settings['control_surface_deflection_generator_settings'][str(i_cs)])
+                        settings['control_surface_deflection_generator_settings'][str(i_cs)])
                 except KeyError:
                     with_error_initialising_cs = True
                     cout.cout_wrap('Error, unable to locate a settings dictionary for control surface '
@@ -135,72 +102,45 @@ def generate(self, aero_dict, beam, aero_settings, ts):
         if with_error_initialising_cs:
             raise KeyError('Unable to locate settings for at least one control surface.')
 
+
         self.add_timestep()
         self.generate_mapping()
         self.generate_zeta(self.beam, self.aero_settings, ts)
 
-        if 'polars' in aero_dict:
+        if 'polars' in self.data_dict:
             import sharpy.aero.utils.airfoilpolars as ap
             self.polars = []
-            nairfoils = np.amax(self.aero_dict['airfoil_distribution']) + 1
+            nairfoils = np.amax(self.data_dict['airfoil_distribution']) + 1
             for iairfoil in range(nairfoils):
                 new_polar = ap.Polar()
-                new_polar.initialise(aero_dict['polars'][str(iairfoil)])
+                new_polar.initialise(data_dict['polars'][str(iairfoil)])
                 self.polars.append(new_polar)
 
     def output_info(self):
         cout.cout_wrap('The aerodynamic grid contains %u surfaces' % self.n_surf, 1)
         for i_surf in range(self.n_surf):
             cout.cout_wrap('  Surface %u, M=%u, N=%u' % (i_surf,
-                                                         self.aero_dimensions[i_surf, 0],
-                                                         self.aero_dimensions[i_surf, 1]), 1)
+                                                         self.dimensions[i_surf, 0],
+                                                         self.dimensions[i_surf, 1]), 1)
             cout.cout_wrap('     Wake %u, M=%u, N=%u' % (i_surf,
-                                                         self.aero_dimensions_star[i_surf, 0],
-                                                         self.aero_dimensions_star[i_surf, 1]))
-        cout.cout_wrap('  In total: %u bound panels' % (sum(self.aero_dimensions[:, 0]*
-                                                            self.aero_dimensions[:, 1])))
-        cout.cout_wrap('  In total: %u wake panels' % (sum(self.aero_dimensions_star[:, 0]*
-                                                           self.aero_dimensions_star[:, 1])))
-        cout.cout_wrap('  Total number of panels = %u' % (sum(self.aero_dimensions[:, 0]*
-                                                              self.aero_dimensions[:, 1]) +
-                                                          sum(self.aero_dimensions_star[:, 0]*
-                                                              self.aero_dimensions_star[:, 1])))
+                                                         self.dimensions_star[i_surf, 0],
+                                                         self.dimensions_star[i_surf, 1]))
+        cout.cout_wrap('  In total: %u bound panels' % (sum(self.dimensions[:, 0]*
+                                                            self.dimensions[:, 1])))
+        cout.cout_wrap('  In total: %u wake panels' % (sum(self.dimensions_star[:, 0]*
+                                                           self.dimensions_star[:, 1])))
+        cout.cout_wrap('  Total number of panels = %u' % (sum(self.dimensions[:, 0]*
+                                                              self.dimensions[:, 1]) +
+                                                          sum(self.dimensions_star[:, 0]*
+                                                              self.dimensions_star[:, 1])))
 
     def calculate_dimensions(self):
-        self.aero_dimensions = np.zeros((self.n_surf, 2), dtype=int)
-        for i in range(self.n_surf):
-            # adding M values
-            self.aero_dimensions[i, 0] = self.surface_m[i]
-        # count N values (actually, the count result
-        # will be N+1)
-        nodes_in_surface = []
-        for i_surf in range(self.n_surf):
-            nodes_in_surface.append([])
-        for i_elem in range(self.beam.num_elem):
-            nodes = self.beam.elements[i_elem].global_connectivities
-            i_surf = self.aero_dict['surface_distribution'][i_elem]
-            if i_surf < 0:
-                continue
-            for i_global_node in nodes:
-                if i_global_node in nodes_in_surface[i_surf]:
-                    continue
-                else:
-                    nodes_in_surface[i_surf].append(i_global_node)
-                if self.aero_dict['aero_node'][i_global_node]:
-                    self.aero_dimensions[i_surf, 1] += 1
-        # accounting for N+1 nodes -> N panels
-        self.aero_dimensions[:,1] -= 1
-        
-        self.aero_dimensions_star = self.aero_dimensions.copy()
-        self.aero_dimensions_star[:, 0] = self.aero_settings['mstar']
+        super().calculate_dimensions()
 
-    def add_timestep(self):
-        try:
-            self.timestep_info.append(self.timestep_info[-1].copy())
-        except IndexError:
-            self.timestep_info.append(self.ini_info.copy())
+        self.dimensions_star = self.dimensions.copy()
+        self.dimensions_star[:, 0] = self.aero_settings['mstar']
 
-    def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_settings, it=None, dt=None):
+    def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, settings, it=None, dt=None):
         if it is None:
             it = len(beam.timestep_info) - 1
         global_node_in_surface = []
@@ -209,24 +149,24 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_se
 
         # check that we have control surface information
         try:
-            self.aero_dict['control_surface']
+            self.data_dict['control_surface']
             with_control_surfaces = True
         except KeyError:
             with_control_surfaces = False
 
         # check that we have sweep information
         try:
-            self.aero_dict['sweep']
+            self.data_dict['sweep']
         except KeyError:
-            self.aero_dict['sweep'] = np.zeros_like(self.aero_dict['twist'])
+            self.data_dict['sweep'] = np.zeros_like(self.data_dict['twist'])
 
         # Define first_twist for backwards compatibility
-        if 'first_twist' not in self.aero_dict:     
-            self.aero_dict['first_twist'] = [True]*self.aero_dict['surface_m'].shape[0]
+        if 'first_twist' not in self.data_dict:     
+            self.data_dict['first_twist'] = [True]*self.data_dict['surface_m'].shape[0]
         
         # one surface per element
         for i_elem in range(self.n_elem):
-            i_surf = self.aero_dict['surface_distribution'][i_elem]
+            i_surf = self.data_dict['surface_distribution'][i_elem]
             # check if we have to generate a surface here
             if i_surf == -1:
                 continue
@@ -235,7 +175,7 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_se
                 i_global_node = self.beam.elements[i_elem].global_connectivities[i_local_node]
                 # i_global_node = self.beam.elements[i_elem].global_connectivities[
                 #     self.beam.elements[i_elem].ordering[i_local_node]]
-                if not self.aero_dict['aero_node'][i_global_node]:
+                if not self.data_dict['aero_node'][i_global_node]:
                     continue
                 if i_global_node in global_node_in_surface[i_surf]:
                     continue
@@ -263,23 +203,23 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_se
                 control_surface_info = None
                 if with_control_surfaces:
                 # 1) check that this node and elem have a control surface
-                    if self.aero_dict['control_surface'][i_elem, i_local_node] >= 0:
-                        i_control_surface = self.aero_dict['control_surface'][i_elem, i_local_node]
-                        # 2) type of control surface + write info
+                    if self.data_dict['control_surface'][i_elem, i_local_node] >= 0:
+                        i_control_surface = self.data_dict['control_surface'][i_elem, i_local_node]
+                # 2) type of control surface + write info
                         control_surface_info = dict()
-                        if self.aero_dict['control_surface_type'][i_control_surface] == 0:
+                        if self.data_dict['control_surface_type'][i_control_surface] == 0:
                             control_surface_info['type'] = 'static'
-                            control_surface_info['deflection'] = self.aero_dict['control_surface_deflection'][i_control_surface]
-                            control_surface_info['chord'] = self.aero_dict['control_surface_chord'][i_control_surface]
+                            control_surface_info['deflection'] = self.data_dict['control_surface_deflection'][i_control_surface]
+                            control_surface_info['chord'] = self.data_dict['control_surface_chord'][i_control_surface]
                             try:
-                                control_surface_info['hinge_coords'] = self.aero_dict['control_surface_hinge_coords'][i_control_surface]
+                                control_surface_info['hinge_coords'] = self.data_dict['control_surface_hinge_coords'][i_control_surface]
                             except KeyError:
                                 control_surface_info['hinge_coords'] = None
-                        elif self.aero_dict['control_surface_type'][i_control_surface] == 1:
+                        elif self.data_dict['control_surface_type'][i_control_surface] == 1:
                             control_surface_info['type'] = 'dynamic'
-                            control_surface_info['chord'] = self.aero_dict['control_surface_chord'][i_control_surface]
+                            control_surface_info['chord'] = self.data_dict['control_surface_chord'][i_control_surface]
                             try:
-                                control_surface_info['hinge_coords'] = self.aero_dict['control_surface_hinge_coords'][i_control_surface]
+                                control_surface_info['hinge_coords'] = self.data_dict['control_surface_hinge_coords'][i_control_surface]
                             except KeyError:
                                 control_surface_info['hinge_coords'] = None
 
@@ -287,7 +227,7 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_se
                             control_surface_info['deflection'], control_surface_info['deflection_dot'] = \
                                 self.cs_generators[i_control_surface](params)
 
-                        elif self.aero_dict['control_surface_type'][i_control_surface] == 2:
+                        elif self.data_dict['control_surface_type'][i_control_surface] == 2:
                             control_surface_info['type'] = 'controlled'
 
                             try:
@@ -296,12 +236,12 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_se
                                 try:
                                     old_deflection = aero_tstep.control_surface_deflection[i_control_surface]
                                 except IndexError:
-                                    old_deflection = self.aero_dict['control_surface_deflection'][i_control_surface]
+                                    old_deflection = self.data_dict['control_surface_deflection'][i_control_surface]
 
                             try:
                                 control_surface_info['deflection'] = aero_tstep.control_surface_deflection[i_control_surface]
                             except IndexError:
-                                control_surface_info['deflection'] = self.aero_dict['control_surface_deflection'][i_control_surface]
+                                control_surface_info['deflection'] = self.data_dict['control_surface_deflection'][i_control_surface]
 
                             if dt is not None:
                                 control_surface_info['deflection_dot'] = (
@@ -309,15 +249,14 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_se
                             else:
                                 control_surface_info['deflection_dot'] = 0.0
 
-                            control_surface_info['chord'] = self.aero_dict['control_surface_chord'][i_control_surface]
+                            control_surface_info['chord'] = self.data_dict['control_surface_chord'][i_control_surface]
 
                             try:
-                                control_surface_info['hinge_coords'] = self.aero_dict['control_surface_hinge_coords'][i_control_surface]
+                                control_surface_info['hinge_coords'] = self.data_dict['control_surface_hinge_coords'][i_control_surface]
                             except KeyError:
                                 control_surface_info['hinge_coords'] = None
-
                         else:
-                            raise NotImplementedError(str(self.aero_dict['control_surface_type'][i_control_surface]) +
+                            raise NotImplementedError(str(self.data_dict['control_surface_type'][i_control_surface]) +
                                 ' control surfaces are not yet implemented')
 
 
@@ -325,13 +264,13 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_se
                 node_info = dict()
                 node_info['i_node'] = i_global_node
                 node_info['i_local_node'] = i_local_node
-                node_info['chord'] = self.aero_dict['chord'][i_elem, i_local_node]
-                node_info['eaxis'] = self.aero_dict['elastic_axis'][i_elem, i_local_node]
-                node_info['twist'] = self.aero_dict['twist'][i_elem, i_local_node]
-                node_info['sweep'] = self.aero_dict['sweep'][i_elem, i_local_node]
-                node_info['M'] = self.aero_dimensions[i_surf, 0]
-                node_info['M_distribution'] = self.aero_dict['m_distribution'].decode('ascii')
-                node_info['airfoil'] = self.aero_dict['airfoil_distribution'][i_elem, i_local_node]
+                node_info['chord'] = self.data_dict['chord'][i_elem, i_local_node]
+                node_info['eaxis'] = self.data_dict['elastic_axis'][i_elem, i_local_node]
+                node_info['twist'] = self.data_dict['twist'][i_elem, i_local_node]
+                node_info['sweep'] = self.data_dict['sweep'][i_elem, i_local_node]
+                node_info['M'] = self.dimensions[i_surf, 0]
+                node_info['M_distribution'] = self.data_dict['m_distribution'].decode('ascii')
+                node_info['airfoil'] = self.data_dict['airfoil_distribution'][i_elem, i_local_node]
                 node_info['control_surface'] = control_surface_info
                 node_info['beam_coord'] = structure_tstep.pos[i_global_node, :]
                 node_info['pos_dot'] = structure_tstep.pos_dot[i_global_node, :]
@@ -343,75 +282,25 @@ def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_se
                 node_info['cga'] = structure_tstep.cga()
                 if node_info['M_distribution'].lower() == 'user_defined':
                     ielem_in_surf = i_elem - np.sum(self.surface_distribution < i_surf)
-                    node_info['user_defined_m_distribution'] = self.aero_dict['user_defined_m_distribution'][str(i_surf)][:, ielem_in_surf, i_local_node]
+                    node_info['user_defined_m_distribution'] = self.data_dict['user_defined_m_distribution'][str(i_surf)][:, ielem_in_surf, i_local_node]
                 (aero_tstep.zeta[i_surf][:, :, i_n],
                  aero_tstep.zeta_dot[i_surf][:, :, i_n]) = (
                     generate_strip(node_info,
                                    self.airfoil_db,
-                                   aero_settings['aligned_grid'],
-                                   orientation_in=aero_settings['freestream_dir'],
-                                   calculate_zeta_dot=True,
-                                   first_twist=self.aero_dict['first_twist'][i_surf]))
-
-    def generate_zeta(self, beam, aero_settings, ts=-1, beam_ts=-1):
-        self.generate_zeta_timestep_info(beam.timestep_info[beam_ts],
-                                         self.timestep_info[ts],
-                                         beam,
-                                         aero_settings)
-
-    def generate_mapping(self):
-        self.struct2aero_mapping = [[]]*self.n_node
-        surf_n_counter = np.zeros((self.n_surf,), dtype=int)
-        nodes_in_surface = []
-        for i_surf in range(self.n_surf):
-            nodes_in_surface.append([])
-
-        for i_elem in range(self.n_elem):
-            i_surf = self.aero_dict['surface_distribution'][i_elem]
-            if i_surf == -1:
-                continue
-            for i_global_node in self.beam.elements[i_elem].reordered_global_connectivities:
-                if not self.aero_dict['aero_node'][i_global_node]:
-                    continue
+                                   self.aero_settings['aligned_grid'],
+                                   orientation_in=self.aero_settings['freestream_dir'],
+                                   calculate_zeta_dot=True))
+        # set junction boundary conditions for later phantom cell creation in UVLM
+        if "junction_boundary_condition" in self.data_dict:
+            if np.any(self.data_dict["junction_boundary_condition"] >= 0):
+                self.generate_phantom_panels_at_junction(aero_tstep)
 
-                if i_global_node in nodes_in_surface[i_surf]:
-                    continue
-                else:
-                    nodes_in_surface[i_surf].append(i_global_node)
-                    surf_n_counter[i_surf] += 1
-                    try:
-                        self.struct2aero_mapping[i_global_node][0]
-                    except IndexError:
-                        self.struct2aero_mapping[i_global_node] = []
-
-                i_n = surf_n_counter[i_surf] - 1
-                self.struct2aero_mapping[i_global_node].append({'i_surf': i_surf,
-                                                                'i_n': i_n})
-
-        nodes_in_surface = []
-        for i_surf in range(self.n_surf):
-            nodes_in_surface.append([])
 
+    def generate_phantom_panels_at_junction(self, aero_tstep):
         for i_surf in range(self.n_surf):
-            self.aero2struct_mapping.append([-1]*(surf_n_counter[i_surf]))
+                aero_tstep.flag_zeta_phantom[0, i_surf] = self.data_dict["junction_boundary_condition"][0,i_surf] 
 
-        for i_elem in range(self.n_elem):
-            for i_global_node in self.beam.elements[i_elem].global_connectivities:
-                for i in range(len(self.struct2aero_mapping[i_global_node])):
-                    try:
-                        i_surf = self.struct2aero_mapping[i_global_node][i]['i_surf']
-                        i_n = self.struct2aero_mapping[i_global_node][i]['i_n']
-                        if i_global_node in nodes_in_surface[i_surf]:
-                            continue
-                        else:
-                            nodes_in_surface[i_surf].append(i_global_node)
-                    except KeyError:
-                        continue
-                    self.aero2struct_mapping[i_surf][i_n] = i_global_node
 
-    def update_orientation(self, quat, ts=-1):
-        rot = algebra.quat2rotation(quat)
-        self.timestep_info[ts].update_orientation(rot.T)
 
     @staticmethod
     def compute_gamma_dot(dt, tstep, previous_tsteps):
@@ -501,6 +390,9 @@ def generate_strip(node_info, airfoil_db, aligned_grid,
         strip_coordinates_b_frame[1, :] = np.linspace(0.0, 1.0, node_info['M'] + 1)
     elif node_info['M_distribution'] == '1-cos':
         domain = np.linspace(0, 1.0, node_info['M'] + 1)
+        strip_coordinates_b_frame[1, :] = 0.5*(1.0 - np.cos(domain*np.pi))
+    elif node_info['M_distribution'].lower() == 'user_defined':
+        strip_coordinates_b_frame[1,:] = node_info['user_defined_m_distribution']
     else:
         raise NotImplemented('M_distribution is ' + node_info['M_distribution'] +
                              ' and it is not yet supported')
@@ -533,7 +425,7 @@ def generate_strip(node_info, airfoil_db, aligned_grid,
             # deflection velocity
             try:
                 cs_velocity[:, i_M] += np.cross(np.array([-node_info['control_surface']['deflection_dot'], 0.0, 0.0]),
-                                                relative_coords)
+                                            relative_coords)
             except KeyError:
                 pass
 
diff --git a/sharpy/aero/models/grid.py b/sharpy/aero/models/grid.py
new file mode 100644
index 000000000..176a7a5e6
--- /dev/null
+++ b/sharpy/aero/models/grid.py
@@ -0,0 +1,174 @@
+"""Grid
+
+Grid contains
+"""
+import ctypes as ct
+import warnings
+
+import numpy as np
+
+import sharpy.utils.algebra as algebra
+import sharpy.utils.generator_interface as gen_interface
+
+
+class Grid(object):
+    """
+    ``Grid``is the parent class for the lifting surface grid and nonlifting
+    body grids.
+
+    It is created by the solver :class:`sharpy.solvers.aerogridloader.AerogridLoader`
+
+    """
+    def __init__(self):
+        self.data_dict = None
+        self.beam = None
+        self.aero_settings = None
+        self.timestep_info = []
+        self.ini_info = None
+
+        self.surface_distribution = None
+        self.surface_m = None
+        self.dimensions = None
+        self.grid_type = None
+
+        self.n_node = 0
+        self.n_elem = 0
+        self.n_surf = 0
+        self.n_aero_node = 0
+        self.grid_type = None
+
+        self.struct2aero_mapping = None
+        self.aero2struct_mapping = []
+
+
+    def generate(self, data_dict, beam, aero_settings, ts):
+
+
+        self.data_dict = data_dict
+        self.beam = beam
+        self.aero_settings = aero_settings
+        # key words = safe in aero_settings? --> grid_type
+        # number of total nodes (structural + aero&struc)
+        self.n_node = len(data_dict[self.grid_type + '_node']) # gridtype + '_node'
+        # number of elements
+        self.n_elem = len(data_dict['surface_distribution'])
+        # surface distribution
+        self.surface_distribution = data_dict['surface_distribution']
+        # number of surfaces
+        temp = set(data_dict['surface_distribution'])
+        #TO-DO: improve: avoid for loops
+        self.n_surf = sum(1 for i in temp if i >= 0)
+        # number of chordwise panels
+        self.surface_m = data_dict['surface_m']
+        # number of aero nodes
+        self.n_aero_node = sum(data_dict[self.grid_type + '_node'])
+
+
+        # get N per surface
+        self.calculate_dimensions()
+
+        # write grid info to screen
+        # self.output_info()
+
+    def calculate_dimensions(self):
+        self.dimensions = np.zeros((self.n_surf, 2), dtype=int)
+        for i in range(self.n_surf):
+            # adding M values
+            self.dimensions[i, 0] = self.surface_m[i]
+        # Improvement:
+        # self.aero.dimensions[:, 0] = self.surface_m[:]
+        # count N values (actually, the count result
+        # will be N+1)
+        nodes_in_surface = []
+
+        #IMPROVEMENT
+        for i_surf in range(self.n_surf):
+            nodes_in_surface.append([])
+
+        # Improvement!
+        for i_elem in range(self.beam.num_elem):
+            nodes = self.beam.elements[i_elem].global_connectivities
+            i_surf = self.surface_distribution[i_elem]
+            if i_surf < 0:
+                continue
+            for i_global_node in nodes:
+                if i_global_node in nodes_in_surface[i_surf]:
+                    continue
+                else:
+                    nodes_in_surface[i_surf].append(i_global_node)
+                if self.data_dict[self.grid_type + '_node'][i_global_node]:
+                    self.dimensions[i_surf, 1] += 1
+
+        # accounting for N+1 nodes -> N panels
+        self.dimensions[:, 1] -= 1
+
+
+    def add_timestep(self):
+        try:
+            self.timestep_info.append(self.timestep_info[-1].copy())
+        except IndexError:
+            self.timestep_info.append(self.ini_info.copy())
+
+    def generate_zeta_timestep_info(self, structure_tstep, aero_tstep, beam, aero_settings, it=None, dt=None):
+        if it is None:
+            it = len(beam.timestep_info) - 1
+
+
+    def generate_zeta(self, beam, aero_settings, ts=-1, beam_ts=-1):
+        self.generate_zeta_timestep_info(beam.timestep_info[beam_ts],
+                                         self.timestep_info[ts],
+                                         beam,
+                                         aero_settings)
+
+    def generate_mapping(self):
+        self.struct2aero_mapping = [[]]*self.n_node
+        surf_n_counter = np.zeros((self.n_surf,), dtype=int)
+        nodes_in_surface = []
+        for i_surf in range(self.n_surf):
+            nodes_in_surface.append([])
+        for i_elem in range(self.n_elem):
+            i_surf = self.surface_distribution[i_elem]
+            if i_surf == -1:
+                continue
+            for i_global_node in self.beam.elements[i_elem].reordered_global_connectivities:
+                if not self.data_dict[self.grid_type + '_node'][i_global_node]:
+                    continue
+
+                if i_global_node in nodes_in_surface[i_surf]:
+                    continue
+                else:
+                    nodes_in_surface[i_surf].append(i_global_node)
+                    surf_n_counter[i_surf] += 1
+                    try:
+                        self.struct2aero_mapping[i_global_node][0]
+                    except IndexError:
+                        self.struct2aero_mapping[i_global_node] = []
+
+                i_n = surf_n_counter[i_surf] - 1
+                self.struct2aero_mapping[i_global_node].append({'i_surf': i_surf,
+                                                                'i_n': i_n})
+
+        nodes_in_surface = []
+        for i_surf in range(self.n_surf):
+            nodes_in_surface.append([])
+
+        for i_surf in range(self.n_surf):
+            self.aero2struct_mapping.append([-1]*(surf_n_counter[i_surf]))
+
+        for i_elem in range(self.n_elem):
+            for i_global_node in self.beam.elements[i_elem].global_connectivities:
+                for i in range(len(self.struct2aero_mapping[i_global_node])):
+                    try:
+                        i_surf = self.struct2aero_mapping[i_global_node][i]['i_surf']
+                        i_n = self.struct2aero_mapping[i_global_node][i]['i_n']
+                        if i_global_node in nodes_in_surface[i_surf]:
+                            continue
+                        else:
+                            nodes_in_surface[i_surf].append(i_global_node)
+                    except KeyError:
+                        continue
+                    self.aero2struct_mapping[i_surf][i_n] = i_global_node
+
+    def update_orientation(self, quat, ts=-1):
+        rot = algebra.quat2rotation(quat)
+        self.timestep_info[ts].update_orientation(rot.T)
diff --git a/sharpy/aero/models/nonliftingbodygrid.py b/sharpy/aero/models/nonliftingbodygrid.py
new file mode 100644
index 000000000..a4c4d9f61
--- /dev/null
+++ b/sharpy/aero/models/nonliftingbodygrid.py
@@ -0,0 +1,121 @@
+"""Nonlifting Body grid
+
+Description
+"""
+
+from sharpy.aero.models.grid import Grid
+from sharpy.utils.datastructures import NonliftingBodyTimeStepInfo
+import numpy as np
+import sharpy.utils.algebra as algebra
+
+
+class NonliftingBodyGrid(Grid):
+    """
+    ``Nonlifting Body Grid`` is the main object containing information of the
+        nonlifting bodygrid, consisting of triangular and quadrilateral panels.
+    It is created by the solver :class:`sharpy.solvers.aerogridloader.AerogridLoader`
+
+    """
+    def __init__(self):
+        super().__init__()
+        self.grid_type = 'nonlifting_body'
+
+
+    def generate(self, data_dict, beam, nonlifting_body_settings, ts):
+        super().generate(data_dict, beam, nonlifting_body_settings, ts)
+
+        # allocating initial grid storage
+        self.ini_info = NonliftingBodyTimeStepInfo(self.dimensions)
+
+        self.add_timestep()
+        self.generate_mapping()
+        self.generate_zeta(self.beam, self.aero_settings, ts)
+
+    def generate_zeta_timestep_info(self, structure_tstep, nonlifting_body_tstep, beam, aero_settings, it=None, dt=None):
+        super().generate_zeta_timestep_info(structure_tstep, nonlifting_body_tstep, beam, aero_settings, it, dt)
+
+        for i_surf in range(self.n_surf):
+            # Get Zeta and Zeta_dot (Panel node positions in G frame)
+            nonlifting_body_tstep.zeta[i_surf], nonlifting_body_tstep.zeta_dot[i_surf] = self.get_zeta_and_zeta_dot(i_surf, structure_tstep)
+
+    def get_zeta_and_zeta_dot(self,i_surf, structure_tstep):
+        numb_radial_nodes = self.dimensions[i_surf][0] +1
+        matrix_nodes = np.zeros((3, numb_radial_nodes,
+                                 self.dimensions[i_surf][1]+1))
+        matrix_nodes_dot = matrix_nodes.copy()
+        array_phi_coordinates = np.linspace(0, 2*np.pi, numb_radial_nodes)
+        # cache sin and cos values
+        array_sin_phi = np.sin(array_phi_coordinates)
+        array_cos_phi = np.cos(array_phi_coordinates)
+
+        cga_rotation_matrix = structure_tstep.cga()
+
+        for node_counter, i_global_node in enumerate(self.aero2struct_mapping[i_surf]):
+            # 1) Set B-Frame position
+            # 1a) get cross-sectional fuselage geometry at node
+            if self.data_dict["shape"].decode() == 'specific':
+                a_ellipse = self.data_dict["a_ellipse"][i_global_node]
+                b_ellipse = self.data_dict["b_ellipse"][i_global_node]
+                z_0 =self.data_dict["z_0_ellipse"][i_global_node]
+                if a_ellipse == 0. or b_ellipse == 0.:
+                    radius = 0
+                else:
+                    radius = a_ellipse*b_ellipse/np.sqrt(
+                        (b_ellipse*array_cos_phi)**2
+                        +(a_ellipse*array_sin_phi)**2)
+            else:
+                radius = self.data_dict["radius"][i_global_node]
+                z_0 = 0
+            
+            # 1b) Get nodes position in B frame
+            matrix_nodes[1, :, node_counter] = radius*array_cos_phi
+            matrix_nodes[2, :, node_counter] = radius*array_sin_phi + z_0
+            
+            # 2) A frame
+            # 2a) Convert structural position from B to A frame 
+            i_elem, i_local_node = self.get_elment_and_local_node_id(i_surf, i_global_node)
+            psi_node = structure_tstep.psi[i_elem, i_local_node,:]
+            if not (psi_node == [0, 0, 0]).all():
+                # just perform roation from B to A if psi not 0
+                Cab = algebra.crv2rotation(psi_node)
+                for idx in range(numb_radial_nodes):
+                    matrix_nodes[:, idx, node_counter] = np.dot(Cab, matrix_nodes[:, idx, node_counter])
+            # 2b) Add beam displacements (expressed in A-frame)
+            for dim in range(3):
+                matrix_nodes[dim, :, node_counter] += structure_tstep.pos[i_global_node,dim]
+
+            # 2c) Add structural beam velocities (expressed in A-frame)
+            for dim in range(3):
+                # velocity due to pos_dot
+                matrix_nodes_dot[dim, :, node_counter] += structure_tstep.pos_dot[i_global_node, dim]
+            # 2d) Add effect of structural beam rotations an node velocity (expressed in A-frame)
+            psi_dot_node = structure_tstep.psi_dot[i_elem, i_local_node,:]
+            omega_a = algebra.crv_dot2omega(psi_node, psi_dot_node)     
+            for idx in range(numb_radial_nodes):
+                matrix_nodes_dot[:, idx, node_counter] += (np.dot(algebra.skew(omega_a), matrix_nodes[:, idx, node_counter]))
+            
+            # 3) Convert position and velocities from A to G frame
+            for idx in range(numb_radial_nodes):
+                matrix_nodes[:, idx, node_counter] = np.dot(cga_rotation_matrix,
+                                          matrix_nodes[:, idx, node_counter])
+                matrix_nodes_dot[:, idx, node_counter] = np.dot(cga_rotation_matrix,
+                                          matrix_nodes_dot[:, idx, node_counter])
+        return matrix_nodes, matrix_nodes_dot
+
+
+    def get_elment_and_local_node_id(self, i_surf, i_global_node):
+        # get beam elements of surface
+        idx_beam_elements_surface = np.where(self.surface_distribution == i_surf)[0]
+        # find element and local node of the global node and return psi
+        for i_elem in idx_beam_elements_surface:
+            if i_global_node in self.beam.elements[i_elem].reordered_global_connectivities:
+                i_local_node = np.where(self.beam.elements[i_elem].reordered_global_connectivities == i_global_node)[0][0]
+                return i_elem, i_local_node
+        raise Exception("The global node %u could not be assigned to any element of surface %u." % (i_global_node, i_surf))
+
+
+
+
+
+
+
diff --git a/sharpy/aero/utils/mapping.py b/sharpy/aero/utils/mapping.py
index 9bde08983..098f088c7 100644
--- a/sharpy/aero/utils/mapping.py
+++ b/sharpy/aero/utils/mapping.py
@@ -11,7 +11,8 @@ def aero2struct_force_mapping(aero_forces,
                               master,
                               conn,
                               cag=np.eye(3),
-                              aero_dict=None):
+                              data_dict=None,
+                              skip_moments_generated_by_forces = False):
     r"""
     Maps the aerodynamic forces at the lattice to the structural nodes
 
@@ -39,7 +40,8 @@ def aero2struct_force_mapping(aero_forces,
         master: Unused
         conn (np.ndarray): Connectivities matrix
         cag (np.ndarray): Transformation matrix between inertial and body-attached reference ``A``
-        aero_dict (dict): Dictionary containing the grid's information.
+        data_dict (dict): Dictionary containing the grid's information.
+        skip_moments_generated_by_forces (bool): Flag to skip local moment calculation.
 
     Returns:
         np.ndarray: structural forces in an ``n_node x 6`` vector
@@ -69,10 +71,18 @@ def aero2struct_force_mapping(aero_forces,
                 cbg = np.dot(cab.T, cag)
 
                 for i_m in range(n_m):
-                    chi_g = zeta[i_surf][:, i_m, i_n] - np.dot(cag.T, pos_def[i_global_node, :])
                     struct_forces[i_global_node, 0:3] += np.dot(cbg, aero_forces[i_surf][0:3, i_m, i_n])
                     struct_forces[i_global_node, 3:6] += np.dot(cbg, aero_forces[i_surf][3:6, i_m, i_n])
-                    struct_forces[i_global_node, 3:6] += np.dot(cbg, algebra.cross3(chi_g, aero_forces[i_surf][0:3, i_m, i_n]))
+                    """
+                        The calculation of the local moment is skipped for fuselage bodies, since this 
+                        leads to noticeably asymmetric aeroforces. This transitional solution makes 
+                        sense since we have only considered the stiff fuselage so far and the pitching 
+                        moment coming from the fuselage is mainly caused by the longitudinal force distribution.
+                        TODO: Correct the calculation of the local moment for the fuselage model.
+                    """
+                    if not skip_moments_generated_by_forces:
+                        chi_g = zeta[i_surf][:, i_m, i_n] - np.dot(cag.T, pos_def[i_global_node, :])
+                        struct_forces[i_global_node, 3:6] += np.dot(cbg, algebra.cross3(chi_g, aero_forces[i_surf][0:3, i_m, i_n]))
 
     return struct_forces
 
diff --git a/sharpy/aero/utils/uvlmlib.py b/sharpy/aero/utils/uvlmlib.py
index 3e0466300..10903b955 100644
--- a/sharpy/aero/utils/uvlmlib.py
+++ b/sharpy/aero/utils/uvlmlib.py
@@ -2,18 +2,16 @@
 import sharpy.utils.ctypes_utils as ct_utils
 
 import ctypes as ct
+from ctypes import *
 import numpy as np
-import platform
-import os
-from sharpy.utils.constants import NDIM, vortex_radius_def
+from sharpy.utils.constants import vortex_radius_def
 
 try:
     UvlmLib = ct_utils.import_ctypes_lib(SharpyDir + '/UVLM', 'libuvlm')
 except OSError:
     UvlmLib = ct_utils.import_ctypes_lib(SharpyDir + '/lib/UVLM/lib', 'libuvlm')
 
-
-
+# TODO: Combine VMOpts and UVMOpts (Class + inheritance)?
 
 class VMopts(ct.Structure):
     """ctypes definition for VMopts class
@@ -26,11 +24,15 @@ class VMopts(ct.Structure):
             bool NewAIC;
             double DelTime;
             bool Rollup;
+            bool only_lifting;
+            bool only_nonlifting;
             unsigned int NumCores;
             unsigned int NumSurfaces;
-            bool cfl1;
+            unsigned int NumSurfacesNonlifting;
             double vortex_radius;
-            double vortex_radius_wake_ind;
+            double vortex_radius_wake_ind;            
+            bool u_ind_by_sources_for_lifting_forces;
+            uint ignore_first_x_nodes_in_force_calculation;
         };
     """
     _fields_ = [("ImageMethod", ct.c_bool),
@@ -40,8 +42,12 @@ class VMopts(ct.Structure):
                 ("NewAIC", ct.c_bool),
                 ("DelTime", ct.c_double),
                 ("Rollup", ct.c_bool),
+                ("only_lifting", ct.c_bool),
+                ("only_nonlifting", ct.c_bool),
+                ("phantom_wing_test", ct.c_bool),
                 ("NumCores", ct.c_uint),
                 ("NumSurfaces", ct.c_uint),
+                ("NumSurfacesNonlifting", ct.c_uint),
                 ("dt", ct.c_double),
                 ("n_rollup", ct.c_uint),
                 ("rollup_tolerance", ct.c_double),
@@ -49,9 +55,12 @@ class VMopts(ct.Structure):
                 ("iterative_solver", ct.c_bool),
                 ("iterative_tol", ct.c_double),
                 ("iterative_precond", ct.c_bool),
-                ("cfl1", ct.c_bool),
                 ("vortex_radius", ct.c_double),
-                ("vortex_radius_wake_ind", ct.c_double)]
+                ("vortex_radius_wake_ind", ct.c_double),
+                ("consider_u_ind_by_sources_for_lifting_forces", ct.c_bool),
+                ("ignore_first_x_nodes_in_force_calculation", ct.c_uint)]
+    
+
 
     def __init__(self):
         ct.Structure.__init__(self)
@@ -62,8 +71,11 @@ def __init__(self):
         self.NewAIC = ct.c_bool(False)  # legacy var
         self.DelTime = ct.c_double(1.0)
         self.Rollup = ct.c_bool(False)
+        self.only_lifting = ct.c_bool(False)
+        self.only_nonlifting = ct.c_bool(False)
         self.NumCores = ct.c_uint(4)
         self.NumSurfaces = ct.c_uint(1)
+        self.NumSurfacesNonlifting = ct.c_uint(0)
         self.dt = ct.c_double(0.01)
         self.n_rollup = ct.c_uint(0)
         self.rollup_tolerance = ct.c_double(1e-5)
@@ -71,26 +83,49 @@ def __init__(self):
         self.iterative_solver = ct.c_bool(False)
         self.iterative_tol = ct.c_double(0)
         self.iterative_precond = ct.c_bool(False)
-        self.cfl1 = ct.c_bool(True)
         self.vortex_radius = ct.c_double(vortex_radius_def)
         self.vortex_radius_wake_ind = ct.c_double(vortex_radius_def)
-        self.rbm_vel_g = np.ctypeslib.as_ctypes(np.zeros((6)))
+        self.phantom_wing_test = ct.c_bool(False)
+        self.consider_u_ind_by_sources_for_lifting_forces = ct.c_bool(False)
+        self.ignore_first_x_nodes_in_force_calculation = ct.c_uint(0)
+
+
+    def set_options(self, options, n_surfaces = 0, n_surfaces_nonlifting = 0):
+        self.Steady = ct.c_bool(True)
+        self.NumSurfaces = ct.c_uint(n_surfaces)
+        self.NumSurfacesNonlifting = ct.c_uint(n_surfaces_nonlifting)
+        self.horseshoe = ct.c_bool(options['horseshoe'])
+        self.dt = ct.c_double(options["rollup_dt"])
+        self.n_rollup = ct.c_uint(options["n_rollup"])
+        self.rollup_tolerance = ct.c_double(options["rollup_tolerance"])
+        self.rollup_aic_refresh = ct.c_uint(options['rollup_aic_refresh'])
+        self.NumCores = ct.c_uint(options['num_cores'])
+        self.iterative_solver = ct.c_bool(options['iterative_solver'])
+        self.iterative_tol = ct.c_double(options['iterative_tol'])
+        self.iterative_precond = ct.c_bool(options['iterative_precond'])
+        self.vortex_radius = ct.c_double(options['vortex_radius'])
+        self.vortex_radius_wake_ind = ct.c_double(options['vortex_radius_wake_ind'])
+
+        self.only_nonlifting = ct.c_bool(options["only_nonlifting"])
+        self.only_lifting = ct.c_bool(not options["nonlifting_body_interactions"])
+        self.phantom_wing_test = ct.c_bool(options["phantom_wing_test"])
+        self.ignore_first_x_nodes_in_force_calculation = ct.c_uint(options["ignore_first_x_nodes_in_force_calculation"])
 
 
 class UVMopts(ct.Structure):
     _fields_ = [("dt", ct.c_double),
                 ("NumCores", ct.c_uint),
                 ("NumSurfaces", ct.c_uint),
-                # ("steady_n_rollup", ct.c_uint),
-                # ("steady_rollup_tolerance", ct.c_double),
-                # ("steady_rollup_aic_refresh", ct.c_uint),
+                ("NumSurfacesNonlifting", ct.c_uint),            
+                ("only_lifting", ct.c_bool),
+                ("only_nonlifting", ct.c_bool),
+                ("phantom_wing_test", ct.c_bool),
                 ("convection_scheme", ct.c_uint),
-                # ("Mstar", ct.c_uint),
                 ("ImageMethod", ct.c_bool),
                 ("iterative_solver", ct.c_bool),
                 ("iterative_tol", ct.c_double),
                 ("iterative_precond", ct.c_bool),
-                ("convect_wake", ct.c_bool),
+                ("convect_wake", ct.c_bool),             
                 ("cfl1", ct.c_bool),
                 ("vortex_radius", ct.c_double),
                 ("vortex_radius_wake_ind", ct.c_double),
@@ -98,15 +133,18 @@ class UVMopts(ct.Structure):
                 ("filter_method", ct.c_uint),
                 ("interp_method", ct.c_uint),
                 ("yaw_slerp", ct.c_double),
-                ("quasi_steady", ct.c_bool),]
+                ("quasi_steady", ct.c_bool),
+                ("num_spanwise_panels_wo_induced_velocity", ct.c_uint),
+                ("consider_u_ind_by_sources_for_lifting_forces", ct.c_bool),
+                ("ignore_first_x_nodes_in_force_calculation", ct.c_uint)]
 
     def __init__(self):
         ct.Structure.__init__(self)
         self.dt = ct.c_double(0.01)
         self.NumCores = ct.c_uint(4)
         self.NumSurfaces = ct.c_uint(1)
+        self.NumSurfacesNonlifting = ct.c_uint(1)
         self.convection_scheme = ct.c_uint(2)
-        # self.Mstar = ct.c_uint(10)
         self.ImageMethod = ct.c_bool(False)
         self.iterative_solver = ct.c_bool(False)
         self.iterative_tol = ct.c_double(0)
@@ -117,6 +155,46 @@ def __init__(self):
         self.vortex_radius_wake_ind = ct.c_double(vortex_radius_def)
         self.yaw_slerp = ct.c_double(0.)
         self.quasi_steady = ct.c_bool(False)
+        self.num_spanwise_panels_wo_induced_velocity = ct.c_uint(0)
+        self.phantom_wing_test = ct.c_bool(False)
+        self.consider_u_ind_by_sources_for_lifting_forces = ct.c_bool(False)
+        self.ignore_first_x_nodes_in_force_calculation = ct.c_uint(0)
+
+    def set_options(self, 
+                    options, 
+                    n_surfaces = 0, 
+                    n_surfaces_nonlifting = 0, 
+                    dt = None, 
+                    convect_wake = False, 
+                    n_span_panels_wo_u_ind = 0,
+                    only_lifting=True):
+        if dt is None:
+            self.dt = ct.c_double(options["dt"])
+        else:
+            self.dt = ct.c_double(dt)
+        self.NumCores = ct.c_uint(options["num_cores"])
+        self.NumSurfaces = ct.c_uint(n_surfaces)
+        self.NumSurfacesNonlifting = ct.c_uint(n_surfaces_nonlifting)
+        self.ImageMethod = ct.c_bool(False)
+        self.convection_scheme = ct.c_uint(options["convection_scheme"])
+        self.iterative_solver = ct.c_bool(options['iterative_solver'])
+        self.iterative_tol = ct.c_double(options['iterative_tol'])
+        self.iterative_precond = ct.c_bool(options['iterative_precond'])
+        self.convect_wake = ct.c_bool(convect_wake)
+        self.cfl1 = ct.c_bool(options['cfl1'])
+        self.vortex_radius = ct.c_double(options['vortex_radius'])
+        self.vortex_radius_wake_ind = ct.c_double(options['vortex_radius_wake_ind'])
+        self.interp_coords = ct.c_uint(options["interp_coords"])
+        self.filter_method = ct.c_uint(options["filter_method"])
+        self.interp_method = ct.c_uint(options["interp_method"])
+        self.yaw_slerp = ct.c_double(options["yaw_slerp"])
+        self.quasi_steady = ct.c_bool(options['quasi_steady'])
+ 
+        self.only_nonlifting = ct.c_bool(options["only_nonlifting"])
+        self.only_lifting = ct.c_bool(only_lifting)
+        self.phantom_wing_test = ct.c_bool(options["phantom_wing_test"])
+        self.ignore_first_x_nodes_in_force_calculation = ct.c_uint(options["ignore_first_x_nodes_in_force_calculation"])
+        self.num_spanwise_panels_wo_induced_velocity = n_span_panels_wo_u_ind
 
 
 class FlightConditions(ct.Structure):
@@ -125,52 +203,29 @@ class FlightConditions(ct.Structure):
                 ("rho", ct.c_double),
                 ("c_ref", ct.c_double)]
 
-    def __init__(self):
+    def __init__(self, rho, vec_u_inf):
         ct.Structure.__init__(self)
+        self.set_flight_conditions(rho, vec_u_inf)
 
-    # def __init__(self, fc_dict):
-    #     ct.Structure.__init__(self)
-    #     self.uinf = fc_dict['FlightCon']['u_inf']
-    #     alpha = fc_dict['FlightCon']['alpha']
-    #     beta = fc_dict['FlightCon']['beta']
-    #     uinf_direction_temp = np.array([1, 0, 0], dtype=ct.c_double)
-    #     self.uinf_direction = np.ctypeslib.as_ctypes(uinf_direction_temp)
-    #     self.rho = fc_dict['FlightCon']['rho_inf']
-    #     self.c_ref = fc_dict['FlightCon']['c_ref']
+    def set_flight_conditions(self, rho, vec_u_inf):
+        self.rho = rho
+        self.uinf = np.ctypeslib.as_ctypes(np.linalg.norm(vec_u_inf))
+        self.uinf_direction = np.ctypeslib.as_ctypes(vec_u_inf/self.uinf)
 
 
 # type for 2d integer matrix
 t_2int = ct.POINTER(ct.c_int)*2
 
-
 def vlm_solver(ts_info, options):
     run_VLM = UvlmLib.run_VLM
-    run_VLM.restype = None
 
     vmopts = VMopts()
-    vmopts.Steady = ct.c_bool(True)
-    vmopts.NumSurfaces = ct.c_uint(ts_info.n_surf)
-    vmopts.horseshoe = ct.c_bool(options['horseshoe'])
-    vmopts.dt = ct.c_double(options["rollup_dt"])
-    vmopts.n_rollup = ct.c_uint(options["n_rollup"])
-    vmopts.rollup_tolerance = ct.c_double(options["rollup_tolerance"])
-    vmopts.rollup_aic_refresh = ct.c_uint(options['rollup_aic_refresh'])
-    vmopts.NumCores = ct.c_uint(options['num_cores'])
-    vmopts.iterative_solver = ct.c_bool(options['iterative_solver'])
-    vmopts.iterative_tol = ct.c_double(options['iterative_tol'])
-    vmopts.iterative_precond = ct.c_bool(options['iterative_precond'])
-    vmopts.cfl1 = ct.c_bool(options['cfl1'])
-    vmopts.vortex_radius = ct.c_double(options['vortex_radius'])
-    vmopts.vortex_radius_wake_ind = ct.c_double(options['vortex_radius_wake_ind'])
-
-    flightconditions = FlightConditions()
-    flightconditions.rho = options['rho']
-    flightconditions.uinf = np.ctypeslib.as_ctypes(np.linalg.norm(ts_info.u_ext[0][:, 0, 0]))
-    flightconditions.uinf_direction = np.ctypeslib.as_ctypes(ts_info.u_ext[0][:, 0, 0]/flightconditions.uinf)
+    vmopts.set_options(options, n_surfaces = ts_info.n_surf)
+
+    flightconditions = FlightConditions(options['rho'], ts_info.u_ext[0][:, 0, 0])
 
     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),
             ct.byref(flightconditions),
@@ -187,97 +242,81 @@ def vlm_solver(ts_info, options):
             p_centre_rot_g)
     ts_info.remove_ctypes_pointers()
 
-
-def uvlm_init(ts_info, options):
-    init_UVLM = UvlmLib.init_UVLM
-    init_UVLM.restype = None
+def vlm_solver_nonlifting_body(ts_info, options):
+    run_linear_source_panel_method = UvlmLib.run_linear_source_panel_method
 
     vmopts = VMopts()
-    vmopts.Steady = ct.c_bool(True)
-    # vmopts.Mstar = ct.c_uint(options['mstar'])
-    vmopts.NumSurfaces = ct.c_uint(ts_info.n_surf)
-    vmopts.horseshoe = ct.c_bool(False)
-    vmopts.dt = options["dt"]
-    try:
-        vmopts.n_rollup = ct.c_uint(options["steady_n_rollup"])
-        vmopts.rollup_tolerance = ct.c_double(options["steady_rollup_tolerance"])
-        vmopts.rollup_aic_refresh = ct.c_uint(options['steady_rollup_aic_refresh'])
-    except KeyError:
-        pass
-    vmopts.NumCores = ct.c_uint(options['num_cores'])
-    vmopts.vortex_radius = ct.c_double(options['vortex_radius'])
-    vmopts.vortex_radius_wake_ind = ct.c_double(options['vortex_radius_wake_ind'])
-    vmopts.quasi_steady = ct.c_bool(options['quasi_steady'])
-
-    flightconditions = FlightConditions()
-    flightconditions.rho = options['rho']
-    flightconditions.uinf = np.ctypeslib.as_ctypes(np.linalg.norm(ts_info.u_ext[0][:, 0, 0]))
-    flightconditions.uinf_direction = np.ctypeslib.as_ctypes(ts_info.u_ext[0][:, 0, 0]/flightconditions.uinf)
-
-    # rbm_vel[0:3] = np.dot(inertial2aero.transpose(), rbm_vel[0:3])
-    # rbm_vel[3:6] = np.dot(inertial2aero.transpose(), rbm_vel[3:6])
-    p_rbm_vel = np.zeros((6,)).ctypes.data_as(ct.POINTER(ct.c_double))
+    vmopts.set_options(options, n_surfaces_nonlifting = ts_info.n_surf)
+
+    flightconditions = FlightConditions(options['rho'], ts_info.u_ext[0][:, 0, 0])
 
     ts_info.generate_ctypes_pointers()
-    init_UVLM(ct.byref(vmopts),
-              ct.byref(flightconditions),
-              ts_info.ct_p_dimensions,
-              ts_info.ct_p_dimensions_star,
-              ts_info.ct_p_u_ext,
-              ts_info.ct_p_zeta,
-              ts_info.ct_p_zeta_star,
-              ts_info.ct_p_zeta_dot,
-              ts_info.ct_p_zeta_star_dot,
-              p_rbm_vel,
-              ts_info.ct_p_gamma,
-              ts_info.ct_p_gamma_star,
-              ts_info.ct_p_normals,
-              ts_info.ct_p_forces)
+    run_linear_source_panel_method(ct.byref(vmopts),
+                       ct.byref(flightconditions),
+                       ts_info.ct_p_dimensions,
+                       ts_info.ct_p_zeta,
+                       ts_info.ct_p_u_ext,
+                       ts_info.ct_p_sigma,
+                       ts_info.ct_p_forces,
+                       ts_info.ct_p_pressure_coefficients)
     ts_info.remove_ctypes_pointers()
 
+def vlm_solver_lifting_and_nonlifting_bodies(ts_info_lifting, ts_info_nonlifting, options):
+    run_VLM_coupled_with_LSPM = UvlmLib.run_VLM_coupled_with_LSPM
+
+    vmopts = VMopts()    
+    vmopts.set_options(options, n_surfaces = ts_info_lifting.n_surf, n_surfaces_nonlifting = ts_info_nonlifting.n_surf)
+
+    flightconditions = FlightConditions(options['rho'], ts_info_lifting.u_ext[0][:, 0, 0])
+
+    p_rbm_vel_g = options['rbm_vel_g'].ctypes.data_as(ct.POINTER(ct.c_double))
+    p_centre_rot = options['centre_rot_g'].ctypes.data_as(ct.POINTER(ct.c_double))
+    ts_info_lifting.generate_ctypes_pointers()
+    ts_info_nonlifting.generate_ctypes_pointers()
+    run_VLM_coupled_with_LSPM(ct.byref(vmopts),
+            ct.byref(flightconditions),
+            ts_info_lifting.ct_p_dimensions,
+            ts_info_lifting.ct_p_dimensions_star,
+            ts_info_lifting.ct_p_zeta,
+            ts_info_lifting.ct_p_zeta_star,
+            ts_info_lifting.ct_p_zeta_dot,
+            ts_info_lifting.ct_p_u_ext,
+            ts_info_lifting.ct_p_gamma,
+            ts_info_lifting.ct_p_gamma_star,
+            ts_info_lifting.ct_p_forces,
+            ts_info_lifting.ct_p_flag_zeta_phantom,
+            ts_info_nonlifting.ct_p_dimensions,
+            ts_info_nonlifting.ct_p_zeta,
+            ts_info_nonlifting.ct_p_u_ext,
+            ts_info_nonlifting.ct_p_sigma,
+            ts_info_nonlifting.ct_p_forces,
+            ts_info_nonlifting.ct_p_pressure_coefficients,
+            p_rbm_vel_g,
+            p_centre_rot)
+
+    ts_info_lifting.remove_ctypes_pointers()
+    ts_info_nonlifting.remove_ctypes_pointers()
+
 
 def uvlm_solver(i_iter, ts_info, struct_ts_info, options, convect_wake=True, dt=None):
+
+    p_rbm_vel = get_ctype_pointer_of_rbm_vel_in_G_frame(struct_ts_info.for_vel.copy(), struct_ts_info.cga())
+    p_centre_rot = options['centre_rot'].ctypes.data_as(ct.POINTER(ct.c_double))
+
     run_UVLM = UvlmLib.run_UVLM
-    run_UVLM.restype = None
 
     uvmopts = UVMopts()
-    if dt is None:
-        uvmopts.dt = ct.c_double(options["dt"])
-    else:
-        uvmopts.dt = ct.c_double(dt)
-    uvmopts.NumCores = ct.c_uint(options["num_cores"])
-    uvmopts.NumSurfaces = ct.c_uint(ts_info.n_surf)
-    uvmopts.ImageMethod = ct.c_bool(False)
-    uvmopts.convection_scheme = ct.c_uint(options["convection_scheme"])
-    uvmopts.iterative_solver = ct.c_bool(options['iterative_solver'])
-    uvmopts.iterative_tol = ct.c_double(options['iterative_tol'])
-    uvmopts.iterative_precond = ct.c_bool(options['iterative_precond'])
-    uvmopts.convect_wake = ct.c_bool(convect_wake)
-    uvmopts.cfl1 = ct.c_bool(options['cfl1'])
-    uvmopts.vortex_radius = ct.c_double(options['vortex_radius'])
-    uvmopts.vortex_radius_wake_ind = ct.c_double(options['vortex_radius_wake_ind'])
-    uvmopts.interp_coords = ct.c_uint(options["interp_coords"])
-    uvmopts.filter_method = ct.c_uint(options["filter_method"])
-    uvmopts.interp_method = ct.c_uint(options["interp_method"])
-    uvmopts.yaw_slerp = ct.c_double(options["yaw_slerp"])
-    uvmopts.quasi_steady = ct.c_bool(options['quasi_steady'])
-
-    flightconditions = FlightConditions()
-    flightconditions.rho = options['rho']
-    flightconditions.uinf = np.ctypeslib.as_ctypes(np.linalg.norm(ts_info.u_ext[0][:, 0, 0]))
-    # direction = np.array([1.0, 0, 0])
-    flightconditions.uinf_direction = np.ctypeslib.as_ctypes(ts_info.u_ext[0][:, 0, 0]/flightconditions.uinf)
-    # flightconditions.uinf_direction = np.ctypeslib.as_ctypes(direction)
-
-    rbm_vel = struct_ts_info.for_vel.copy()
-    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))
+    uvmopts.set_options(options,
+                        n_surfaces = ts_info.n_surf,
+                        n_surfaces_nonlifting = 0, 
+                        dt = dt, 
+                        convect_wake = convect_wake, 
+                        n_span_panels_wo_u_ind=0)
+
+    flightconditions = FlightConditions(options['rho'], ts_info.u_ext[0][:, 0, 0])
 
     i = ct.c_uint(i_iter)
     ts_info.generate_ctypes_pointers()
-    # previous_ts_info.generate_ctypes_pointers()
     run_UVLM(ct.byref(uvmopts),
              ct.byref(flightconditions),
              ts_info.ct_p_dimensions,
@@ -288,18 +327,64 @@ def uvlm_solver(i_iter, ts_info, struct_ts_info, options, convect_wake=True, dt=
              ts_info.ct_p_zeta,
              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,
-             # previous_ts_info.ct_p_gamma,
              ts_info.ct_p_normals,
              ts_info.ct_p_forces,
-             ts_info.ct_p_dynamic_forces)
+             ts_info.ct_p_dynamic_forces,
+             p_rbm_vel,
+             p_centre_rot)
     ts_info.remove_ctypes_pointers()
-    # previous_ts_info.remove_ctypes_pointers()
 
+def uvlm_solver_lifting_and_nonlifting(i_iter, ts_info, ts_info_nonlifting, struct_ts_info, options, convect_wake=True, dt=None):
+
+    p_rbm_vel = get_ctype_pointer_of_rbm_vel_in_G_frame(struct_ts_info.for_vel.copy(), struct_ts_info.cga())
+    p_centre_rot = options['centre_rot'].ctypes.data_as(ct.POINTER(ct.c_double))
+
+    uvmopts = UVMopts()
+    uvmopts.set_options(options,
+                        n_surfaces = ts_info.n_surf,
+                        n_surfaces_nonlifting = ts_info_nonlifting.n_surf, 
+                        dt = dt, 
+                        convect_wake = convect_wake, 
+                        n_span_panels_wo_u_ind=4,
+                        only_lifting=False)
+    run_UVLM = UvlmLib.run_UVLM_coupled_with_LSPM
+
+    flightconditions = FlightConditions(options['rho'], ts_info.u_ext[0][:, 0, 0])
+
+    i = ct.c_uint(i_iter)
+    ts_info.generate_ctypes_pointers()
+    ts_info_nonlifting.generate_ctypes_pointers()
+    run_UVLM(ct.byref(uvmopts),
+             ct.byref(flightconditions),
+             ts_info.ct_p_dimensions,
+             ts_info.ct_p_dimensions_star,
+             ct.byref(i),
+             ts_info.ct_p_u_ext,
+             ts_info.ct_p_u_ext_star,
+             ts_info.ct_p_zeta,
+             ts_info.ct_p_zeta_star,
+             ts_info.ct_p_zeta_dot,
+             ts_info.ct_p_gamma,
+             ts_info.ct_p_gamma_star,
+             ts_info.ct_p_dist_to_orig,
+             ts_info.ct_p_normals,
+             ts_info.ct_p_forces,
+             ts_info.ct_p_dynamic_forces,
+             ts_info.ct_p_flag_zeta_phantom,
+             ts_info_nonlifting.ct_p_dimensions,
+             ts_info_nonlifting.ct_p_zeta,
+             ts_info_nonlifting.ct_p_u_ext,
+             ts_info_nonlifting.ct_p_sigma,
+             ts_info_nonlifting.ct_p_forces,
+             ts_info_nonlifting.ct_p_pressure_coefficients,
+             p_rbm_vel,
+             p_centre_rot)
+
+    ts_info.remove_ctypes_pointers()
+    ts_info_nonlifting.remove_ctypes_pointers()
 
 def uvlm_calculate_unsteady_forces(ts_info,
                                    struct_ts_info,
@@ -307,7 +392,6 @@ def uvlm_calculate_unsteady_forces(ts_info,
                                    convect_wake=True,
                                    dt=None):
     calculate_unsteady_forces = UvlmLib.calculate_unsteady_forces
-    calculate_unsteady_forces.restype = None
 
     uvmopts = UVMopts()
     if dt is None:
@@ -324,15 +408,10 @@ def uvlm_calculate_unsteady_forces(ts_info,
     uvmopts.convect_wake = ct.c_bool(convect_wake)
     uvmopts.vortex_radius = ct.c_double(options['vortex_radius'])
 
-    flightconditions = FlightConditions()
-    flightconditions.rho = options['rho']
-    flightconditions.uinf = np.ctypeslib.as_ctypes(np.linalg.norm(ts_info.u_ext[0][:, 0, 0]))
-    flightconditions.uinf_direction = np.ctypeslib.as_ctypes(ts_info.u_ext[0][:, 0, 0]/flightconditions.uinf)
 
-    rbm_vel = struct_ts_info.for_vel.copy()
-    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))
+    flightconditions = FlightConditions(options['rho'], ts_info.u_ext[0][:, 0, 0])
+
+    p_rbm_vel = get_ctype_pointer_of_rbm_vel_in_G_frame(struct_ts_info.for_vel.copy(), struct_ts_info.cga())
 
     for i_surf in range(ts_info.n_surf):
         ts_info.dynamic_forces[i_surf].fill(0.0)
@@ -356,12 +435,8 @@ def uvlm_calculate_unsteady_forces(ts_info,
 def uvlm_calculate_incidence_angle(ts_info,
                                    struct_ts_info):
     calculate_incidence_angle = UvlmLib.UVLM_check_incidence_angle
-    calculate_incidence_angle.restype = None
 
-    rbm_vel = struct_ts_info.for_vel.copy()
-    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_rbm_vel = get_ctype_pointer_of_rbm_vel_in_G_frame(struct_ts_info.for_vel.copy(), struct_ts_info.cga())
 
     n_surf = ct.c_uint(ts_info.n_surf)
 
@@ -399,7 +474,6 @@ def uvlm_calculate_total_induced_velocity_at_points(ts_info,
 
     """
     calculate_uind_at_points = UvlmLib.total_induced_velocity_at_points
-    calculate_uind_at_points.restype = None
 
     uvmopts = UVMopts()
     uvmopts.NumSurfaces = ct.c_uint(ts_info.n_surf)
@@ -548,7 +622,6 @@ def get_induced_velocity_cpp(maps, zeta, gamma, zeta_target,
 
     """
     call_ind_vel = UvlmLib.call_ind_vel
-    call_ind_vel.restype = None
 
     assert zeta_target.flags['C_CONTIGUOUS'], "Input not C contiguous"
 
@@ -663,3 +736,9 @@ def dvinddzeta_cpp(zetac, surf_in, is_bound,
     )
 
     return der_coll, der_vert
+
+def get_ctype_pointer_of_rbm_vel_in_G_frame(rbm_vel, cga):
+    rbm_vel[0:3] = np.dot(cga, rbm_vel[0:3])
+    rbm_vel[3:6] = np.dot(cga, rbm_vel[3:6])
+    p_rbm_vel = rbm_vel.ctypes.data_as(ct.POINTER(ct.c_double))
+    return p_rbm_vel
\ No newline at end of file
diff --git a/sharpy/cases/templates/fuselage_wing_configuration/__init__.py b/sharpy/cases/templates/fuselage_wing_configuration/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/sharpy/cases/templates/fuselage_wing_configuration/fuselage_wing_configuration.py b/sharpy/cases/templates/fuselage_wing_configuration/fuselage_wing_configuration.py
new file mode 100644
index 000000000..9f0db4440
--- /dev/null
+++ b/sharpy/cases/templates/fuselage_wing_configuration/fuselage_wing_configuration.py
@@ -0,0 +1,69 @@
+import configobj
+from .fwc_structure import FWC_Structure
+from .fwc_aero import FWC_Aero
+from .fwc_fuselage import FWC_Fuselage
+import os
+import sharpy.sharpy_main
+
+
+class Fuselage_Wing_Configuration:
+    """
+        Fuselage_Wing_Configuration is a template class to create an 
+        aircraft model of a simple wing-fuselage configuration within
+        SHARPy. 
+       
+    """
+
+    def __init__(self, case_name, case_route, output_route):
+        self.case_name = case_name
+        self.case_route = case_route
+        self.output_route = output_route
+
+        self.structure = None
+        self.aero = None
+        self.fuselage = None
+
+        self.settings = None
+
+    def init_aeroelastic(self, **kwargs):
+        self.clean()
+        self.init_structure(**kwargs)
+        self.init_aero(**kwargs)
+        if not kwargs.get('lifting_only', True):
+            self.init_fuselage(**kwargs)
+
+    def init_structure(self, **kwargs):
+        self.structure = FWC_Structure(self.case_name, self.case_route, **kwargs)
+
+    def init_aero(self, **kwargs):
+        self.aero = FWC_Aero(self.structure, self.case_name, self.case_route,**kwargs)
+
+    def init_fuselage(self, **kwargs):
+        self.fuselage = FWC_Fuselage(self.structure, self.case_name, self.case_route,**kwargs)
+    def generate(self):
+        if not os.path.isdir(self.case_route):
+            os.makedirs(self.case_route)
+        self.structure.generate()
+        self.aero.generate()
+        if self.fuselage is not None:
+            self.fuselage.generate()
+
+    def create_settings(self, settings):
+        file_name = self.case_route + '/' + self.case_name + '.sharpy'
+        config = configobj.ConfigObj()
+        config.filename = file_name
+        for k, v in settings.items():
+            config[k] = v
+        config.write()
+        self.settings = settings
+
+    def clean(self):
+        list_files = ['.fem.h5', '.aero.h5', '.nonlifting_body.h5', '.dyn.h5', '.mb.h5', '.sharpy', '.flightcon.txt']
+        for file in list_files:
+            path_file = self.case_route + '/' + self.case_name + file
+            if os.path.isfile(path_file):
+                os.remove(path_file)
+
+    def run(self):
+        sharpy.sharpy_main.main(['', self.case_route + '/' + self.case_name + '.sharpy'])
+
diff --git a/sharpy/cases/templates/fuselage_wing_configuration/fwc_aero.py b/sharpy/cases/templates/fuselage_wing_configuration/fwc_aero.py
new file mode 100644
index 000000000..685f34b71
--- /dev/null
+++ b/sharpy/cases/templates/fuselage_wing_configuration/fwc_aero.py
@@ -0,0 +1,134 @@
+import h5py as h5
+import numpy as np
+
+
+class FWC_Aero:
+    """
+        FWC_Aero contains all attributes to define the aerodynamic grid and makes it accessible for SHARPy.
+    """
+    def __init__(self, structure, case_name, case_route, **kwargs):
+        """        
+        Key-Word Arguments:
+            - structure: structural object of BFF model
+        """
+        self.structure = structure
+
+        self.route = case_route
+        self.case_name = case_name
+
+        self.ea_wing = kwargs.get('elastic_axis',0.5)
+        self.num_chordwise_panels = kwargs.get('num_chordwise_panels', 4)
+        self.chord_wing = kwargs.get('chord', 1.)
+        self.alpha_zero_deg = kwargs.get('alpha_zero_deg', 0.)
+        self.n_surfaces = 2
+        self.radius_fuselage = kwargs.get('max_radius', 0.5)
+        self.lifting_only = kwargs.get('lifting_only', True)
+        
+
+    def generate(self):
+        """
+            Function to set up all necessary parameter inputs to define the geometry and discretisation
+            of the lifting surfaces. Finally, informations are saved to an .h5 file serving as an input
+            file for SHARPy.
+        """
+        self.initialize_attributes()
+        self.set_wing_properties()
+        self.set_junction_boundary_conditions()
+        self.write_input_file()
+
+    def initialize_attributes(self):
+        """
+            Initilializes all necessary attributes for the aero.h5 input file based on the number of 
+            nodes, elements, and surfaces of the model.
+        """
+        self.airfoil_distribution = np.zeros((self.structure.n_elem, self.structure.n_node_elem), dtype=int)
+        self.surface_distribution = np.zeros((self.structure.n_elem,), dtype=int)
+        self.surface_m = np.zeros((self.n_surfaces, ), dtype=int) + self.num_chordwise_panels
+        self.aero_node = np.zeros((self.structure.n_node,), dtype=bool)
+
+        self.twist = np.zeros((self.structure.n_elem, self.structure.n_node_elem))
+        self.sweep = np.zeros_like(self.twist)
+        self.chord = np.zeros_like(self.twist)
+        self.elastic_axis = np.zeros_like(self.twist)
+
+        self.junction_boundary_condition_aero = np.zeros((1, self.n_surfaces), dtype=int) - 1
+    
+    def get_y_junction(self):
+        if self.structure.vertical_wing_position == 0:
+            return self.radius_fuselage
+        else:
+            return np.sqrt(self.radius_fuselage**2-self.structure.vertical_wing_position**2)
+    def set_wing_properties(self):
+        """
+            Sets necessary parameters to define the lifting surfaces of one wing (right).
+        """
+        
+        if not self.lifting_only:
+            self.aero_node[:self.structure.n_node_right_wing] = abs(self.structure.y[:self.structure.n_node_right_wing]) > self.get_y_junction() - 0.05
+            self.aero_node[self.structure.n_node_right_wing:self.structure.n_node_wing_total] = self.aero_node[1:self.structure.n_node_right_wing]
+        else:
+            self.aero_node[:self.structure.n_node_wing_total] = True
+        self.chord[:2*self.structure.n_elem_per_wing, :] = self.chord_wing
+        self.elastic_axis[:2*self.structure.n_elem_per_wing, :] = self.ea_wing
+        self.twist[:2*self.structure.n_elem_per_wing, :] = -np.deg2rad(self.alpha_zero_deg)
+        # surf distribution 0 for right and 1 for left wing
+        self.surface_distribution[self.structure.n_elem_per_wing:2*self.structure.n_elem_per_wing] = 1
+
+    def set_junction_boundary_conditions(self):
+        """
+            Sets the boundary conditions for the fuselage-wing junction. These BCs
+            define the partner surface.
+        """
+        # Right wing (surface 0) has the left wing (surface 1) as a partner surface.
+        self.junction_boundary_condition_aero[0, 0] = 1
+        # Left wing (surface 1) has the left wing (surface 0) as a partner surface.
+        self.junction_boundary_condition_aero[0, 1] = 0
+
+    def write_input_file(self):
+        """
+            Writes previously defined parameters to an .h5 file which serves later as an 
+            input file for SHARPy.
+        """
+            
+        with h5.File(self.route + '/' + self.case_name + '.aero.h5', 'a') as h5file:
+            airfoils_group = h5file.create_group('airfoils')
+            # add one airfoil
+            airfoils_group.create_dataset('0', 
+                                          data=np.column_stack(
+                                            self.generate_naca_camber(P=0, M=0)
+                                          ))
+
+            h5file.create_dataset('chord', data=self.chord)
+            h5file.create_dataset('twist', data=self.twist)
+            h5file.create_dataset('sweep', data=self.sweep)
+
+            # airfoil distribution
+            h5file.create_dataset('airfoil_distribution', data=self.airfoil_distribution)
+            h5file.create_dataset('surface_distribution', data=self.surface_distribution)
+            h5file.create_dataset('surface_m', data=self.surface_m)
+            h5file.create_dataset('m_distribution', data='uniform')
+            h5file.create_dataset('aero_node', data=self.aero_node)
+            h5file.create_dataset('elastic_axis', data=self.elastic_axis)
+            h5file.create_dataset('junction_boundary_condition', data=self.junction_boundary_condition_aero)
+   
+
+    def generate_naca_camber(self,M=0, P=0):
+        """
+            Generates the camber line coordinates of a specified NACA airfoil.
+            # TODO: needed?
+        """
+        mm = M*1e-2
+        p = P*1e-1
+
+        def naca(x, mm, p):
+            if x < 1e-6:
+                return 0.0
+            elif x < p:
+                return mm/(p*p)*(2*p*x - x*x)
+            elif x > p and x < 1+1e-6:
+                return mm/((1-p)*(1-p))*(1 - 2*p + 2*p*x - x*x)
+
+        x_vec = np.linspace(0, 1, 1000)
+        y_vec = np.array([naca(x, mm, p) for x in x_vec])
+        
+        return x_vec, y_vec
diff --git a/sharpy/cases/templates/fuselage_wing_configuration/fwc_fuselage.py b/sharpy/cases/templates/fuselage_wing_configuration/fwc_fuselage.py
new file mode 100644
index 000000000..fd0d30a08
--- /dev/null
+++ b/sharpy/cases/templates/fuselage_wing_configuration/fwc_fuselage.py
@@ -0,0 +1,126 @@
+#! /usr/bin/env python3
+import h5py as h5
+import numpy as np
+import matplotlib.pyplot as plt
+   
+
+class FWC_Fuselage:
+    """
+        FWC_Fuselage contains all attributes to define the nonlifting body grid 
+        representing the fuselage shape and makes it accessible for SHARPy.
+    """
+    def __init__(self, structure, case_name, case_route, **kwargs):
+        """        
+        Key-Word Arguments:
+            - structure: structural object of BFF model
+        """
+        self.num_radial_panels = kwargs.get('num_radial_panels', 12)
+        self.max_radius = kwargs.get('max_radius')
+        self.fuselage_shape = kwargs.get('fuselage_shape','cylindrical')
+        self.flag_plot_radius = kwargs.get('flag_plot_radius', False)
+        self.structure = structure
+
+        self.route = case_route
+        self.case_name = case_name
+        self.n_nonlifting_bodies = 1
+
+    def generate(self): 
+        """
+            Function to set up all necessary parameter inputs to define the geometry and discretisation
+            of the nonlifting surfaces. Finally, informations are saved to an .h5 file serving as an input
+            file for SHARPy.
+        """
+
+        self.initialize_parameters()
+
+        self.nonlifting_body_node[self.structure.n_node_wing_total:self.structure.n_node_fuselage_tail] = True
+        if self.structure.vertical_wing_position == 0:
+            self.nonlifting_body_node[0] = True
+        self.nonlifting_body_distribution[self.structure.n_elem_per_wing*2:self.structure.n_elem_per_wing*2+self.structure.n_elem_fuselage] = 0
+        self.nonlifting_body_m[0] = self.num_radial_panels
+
+        self.get_fuselage_geometry()
+
+        self.write_fuselage_input_file()
+        
+
+    def initialize_parameters(self):  
+        """
+            Initilializes all necessary attributes for the nonlifting.h5 input file based on the number of 
+            nodes, elements, and surfaces of the nonlifting model.
+        """      
+        self.nonlifting_body_node = np.zeros((self.structure.n_node,), dtype=bool)
+        self.nonlifting_body_distribution = np.zeros((self.structure.n_elem,), dtype=int) - 1
+        self.nonlifting_body_m = np.zeros((self.n_nonlifting_bodies, ), dtype=int)
+        self.radius = np.zeros((self.structure.n_node,))
+
+    def get_fuselage_geometry(self):
+        x_coord_fuselage_sorted = np.sort(self.get_values_at_fuselage_nodes(self.structure.x))
+        if self.fuselage_shape == 'cylindrical':
+            radius_fuselage = self.create_fuselage_geometry(x_coord_fuselage_sorted.copy(), 0.2*self.structure.fuselage_length, 0.8*self.structure.fuselage_length) 
+        elif self.fuselage_shape == 'ellipsoid':
+            radius_fuselage = self.get_radius_ellipsoid(x_coord_fuselage_sorted.copy(), self.structure.fuselage_length/2, self.max_radius)
+        else:
+            raise "ERROR Fuselage shape {} unknown.".format(self.fuselage_shape)
+        if self.structure.vertical_wing_position == 0:
+            self.radius[0] = radius_fuselage[self.structure.idx_junction]
+            self.radius[self.structure.n_node_wing_total:] =np.delete(radius_fuselage, self.structure.idx_junction)
+        else:
+            self.radius[self.structure.n_node_wing_total:self.structure.n_node_fuselage_tail] = radius_fuselage
+        if self.flag_plot_radius:
+            self.plot_fuselage_radius(self.get_values_at_fuselage_nodes(self.structure.x), self.get_values_at_fuselage_nodes(self.radius))
+            self.plot_fuselage_radius(x_coord_fuselage_sorted, radius_fuselage)
+
+    def get_values_at_fuselage_nodes(self, array):
+        return array[self.nonlifting_body_node]
+    
+    def plot_fuselage_radius(self, x, radius):
+        plt.scatter(x,
+                 radius)
+        plt.grid()
+        plt.xlabel("x, m")
+        plt.ylabel("r, m")
+        plt.gca().set_aspect('equal')
+        plt.show()
+
+    def write_fuselage_input_file(self):
+        with h5.File(self.route + '/' + self.case_name + '.nonlifting_body.h5', 'a') as h5file:
+            h5file.create_dataset('shape', data='cylindrical')
+            h5file.create_dataset('surface_m', data=self.nonlifting_body_m)
+            h5file.create_dataset('nonlifting_body_node', data=self.nonlifting_body_node)
+            h5file.create_dataset('surface_distribution', data=self.nonlifting_body_distribution)
+            h5file.create_dataset('radius', data=self.radius)
+
+    def get_radius_ellipsoid(self, x_coordinates, a, b):
+        x_coordinates[:] -=  x_coordinates[-1] - a
+        y_coordinates = b*np.sqrt(1-(x_coordinates/a)**2)
+        return y_coordinates
+
+    def find_index_of_closest_entry(self, array_values, target_value):
+        return np.argmin(np.abs(array_values - target_value))
+
+
+    def create_fuselage_geometry(self, x_coord_fuselage, x_nose_end, x_tail_start):
+        array_radius = np.zeros_like(x_coord_fuselage)
+        # start with nose at zero to get indices of nose and tail correct
+        x_coord_fuselage[:] -=  x_coord_fuselage[0]
+        idx_cylinder_start = self.find_index_of_closest_entry(x_coord_fuselage, x_nose_end)
+        idx_cylinder_end = self.find_index_of_closest_entry(x_coord_fuselage, x_tail_start)
+
+        # set constant radius of cylinder
+        array_radius[idx_cylinder_start:idx_cylinder_end] = self.max_radius
+        # get ellipsoidal nose and tail shape
+        array_radius[:idx_cylinder_start+1] = self.get_nose_shape(x_coord_fuselage[:idx_cylinder_start+1],
+                                                                  self.max_radius)
+        
+        array_radius[idx_cylinder_end-1:] = np.flip(self.get_nose_shape(x_coord_fuselage[idx_cylinder_end-1:],
+                                                                  self.max_radius))
+        
+        return array_radius
+    
+    def get_nose_shape(self, x_coord, radius):
+        n_nodes = len(x_coord)
+        x_coord[:] -= x_coord[-1]
+        x_coord = np.concatenate((x_coord, np.flip(x_coord[:-1])))        
+        radius = self.get_radius_ellipsoid(x_coord, x_coord[0], radius)
+        return radius[:n_nodes]
diff --git a/sharpy/cases/templates/fuselage_wing_configuration/fwc_get_settings.py b/sharpy/cases/templates/fuselage_wing_configuration/fwc_get_settings.py
new file mode 100644
index 000000000..311ff1ac6
--- /dev/null
+++ b/sharpy/cases/templates/fuselage_wing_configuration/fwc_get_settings.py
@@ -0,0 +1,153 @@
+import numpy as np
+import sharpy.utils.algebra as algebra
+
+def define_simulation_settings(flow, model, alpha_deg, u_inf, 
+                               dt=1,
+                               rho = 1.225, 
+                               lifting_only=True, 
+                               nonlifting_only=False, 
+                               phantom_test=False,
+                               horseshoe=False, **kwargs):
+    gravity = kwargs.get('gravity',True)
+    nonlifting_body_interactions = not lifting_only and not nonlifting_only
+    wake_length = kwargs.get('wake_length', 10)
+    # Other parameters
+    if horseshoe:
+        mstar = 1 
+    else:
+        mstar = wake_length*model.aero.num_chordwise_panels
+    # numerics
+    n_step = kwargs.get('n_step', 5)
+    structural_relaxation_factor = kwargs.get('structural_relaxation_factor', 0.6)
+    tolerance = kwargs.get('tolerance', 1e-6)
+    fsi_tolerance = kwargs.get('fsi_tolerance', 1e-6)
+    num_cores = kwargs.get('num_cores',2)
+
+    if not lifting_only:
+        nonlifting_body_interactions = True
+    else:
+        nonlifting_body_interactions = False
+    settings = {}
+    settings['SHARPy'] = {'case': model.case_name,
+                        'route': model.case_route,
+                        'flow': flow,
+                        'write_screen': 'off',
+                        'write_log': 'on',
+                        'log_folder': model.output_route,
+                        'log_file': model.case_name + '.log'}
+
+
+    settings['BeamLoader'] = {'unsteady': 'on',
+                                'orientation': algebra.euler2quat(np.array([0.,
+                                                                            np.deg2rad(alpha_deg),
+                                                                            0.]))}
+
+    settings['BeamLoads']  = {}
+    settings['LiftDistribution'] = {'coefficients': True}
+    
+    settings['NonLinearStatic'] = {'print_info': 'off',
+                                'max_iterations': 150,
+                                'num_load_steps': 1,
+                                'delta_curved': 1e-1,
+                                'min_delta': tolerance,
+                                'gravity_on': gravity,
+                                'gravity': 9.81}
+
+    settings['StaticUvlm'] = {'print_info': 'on',
+                            'horseshoe': horseshoe,
+                            'num_cores': num_cores,
+                            'velocity_field_generator': 'SteadyVelocityField',
+                            'velocity_field_input': {'u_inf': u_inf,
+                                                    'u_inf_direction': [1., 0, 0]},
+                            'rho': rho,
+                            'nonlifting_body_interactions': nonlifting_body_interactions,
+                            'only_nonlifting': nonlifting_only,
+                            'phantom_wing_test': phantom_test,
+                            'ignore_first_x_nodes_in_force_calculation': kwargs.get('ignore_first_x_nodes_in_force_calculation', 0),
+                            }
+
+    settings['StaticCoupled'] = {'print_info': 'off',
+                                'structural_solver': 'NonLinearStatic',
+                                'structural_solver_settings': settings['NonLinearStatic'],
+                                'aero_solver': 'StaticUvlm',
+                                'aero_solver_settings': settings['StaticUvlm'],
+                                'max_iter': 100,
+                                'n_load_steps': n_step,
+                                'tolerance': fsi_tolerance,
+                                'relaxation_factor': structural_relaxation_factor,
+                                'nonlifting_body_interactions': nonlifting_body_interactions}
+
+    settings['AerogridLoader'] = {'unsteady': 'on',
+                                'aligned_grid': 'on',
+                                'mstar': mstar, #int(20/tstep_factor),
+                                'wake_shape_generator': 'StraightWake',
+                                'wake_shape_generator_input': {
+                                    'u_inf': u_inf,
+                                    'u_inf_direction': [1., 0., 0.],
+                                    'dt': dt,
+                                },
+                            }
+    if 'WriteVariablesTime' in flow:
+        settings['WriteVariablesTime'] = {'cleanup_old_solution': True}
+        if kwargs.get('writeCpVariables', False):
+            settings['WriteVariablesTime']['nonlifting_nodes_variables'] = ['pressure_coefficients']                                                           
+            settings['WriteVariablesTime']['nonlifting_nodes_isurf'] = np.zeros((model.structure.n_node_fuselage,))
+            settings['WriteVariablesTime']['nonlifting_nodes_im'] =  np.zeros((model.structure.n_node_fuselage))
+            settings['WriteVariablesTime']['nonlifting_nodes_in'] = list(range(model.structure.n_node_fuselage))
+        if kwargs.get('writeWingPosVariables', False):
+            settings['WriteVariablesTime']['structure_variables'] = ['pos']
+            settings['WriteVariablesTime']['structure_nodes'] = list(range(model.structure.n_node-2))
+
+    settings['NonliftingbodygridLoader'] = {}
+
+    settings['AeroForcesCalculator'] = {'coefficients': False}
+    settings['AerogridPlot'] = {'plot_nonlifting_surfaces': nonlifting_body_interactions}
+    settings['BeamPlot'] = {}
+
+    settings['NoStructural'] = {}
+    settings['NonLinearDynamicPrescribedStep'] = {'print_info': 'off',
+                                            'max_iterations': 950,
+                                            'delta_curved': 1e-1,
+                                            'min_delta': 1e-6,
+                                            'newmark_damp': 1e-4,
+                                            'gravity_on': gravity,
+                                            'gravity': 9.81,
+                                            'num_steps': kwargs.get('n_tsteps',10),
+                                            'dt': dt,
+                                            }
+    settings['StepUvlm'] = {'print_info': 'on',
+                                'num_cores': 4,
+                                'convection_scheme': 3,
+                                'velocity_field_input': {'u_inf': u_inf,
+                                                            'u_inf_direction': [1., 0., 0.]},
+                                'rho': rho,
+                                'n_time_steps': kwargs.get('n_tsteps',10),
+                                'dt': dt,
+                                'phantom_wing_test': phantom_test,
+                                'nonlifting_body_interactions': not lifting_only,
+                                'gamma_dot_filtering': 3,                                
+                                'ignore_first_x_nodes_in_force_calculation': kwargs.get('ignore_first_x_nodes_in_force_calculation', 0),}
+    
+    dynamic_structural_solver = kwargs.get('structural_solver','NonLinearDynamicPrescribedStep')
+    settings['DynamicCoupled'] = {'structural_solver': dynamic_structural_solver,
+                                    'structural_solver_settings': settings[dynamic_structural_solver],
+                                    'aero_solver': 'StepUvlm',
+                                    'aero_solver_settings': settings['StepUvlm'],
+                                    'fsi_substeps': kwargs.get('fsi_substeps', 200),
+                                    'fsi_tolerance': fsi_tolerance,
+                                    'relaxation_factor': kwargs.get('relaxation_factor',0.1),
+                                    'minimum_steps': 1,
+                                    'relaxation_steps': 150,
+                                    'final_relaxation_factor': kwargs.get('final_relaxation_factor', 0.05),
+                                    'n_time_steps': kwargs.get('n_tsteps',10),
+                                    'dt': dt,
+                                    'nonlifting_body_interactions': not lifting_only,
+                                    'include_unsteady_force_contribution': kwargs.get('unsteady_force_distribution', True), 
+                                    'postprocessors': ['BeamLoads'],
+                                    'postprocessors_settings': {
+                                                                'BeamLoads': {'csv_output': 'off'},
+                                                                },
+                                }
+
+
+    return settings
\ No newline at end of file
diff --git a/sharpy/cases/templates/fuselage_wing_configuration/fwc_structure.py b/sharpy/cases/templates/fuselage_wing_configuration/fwc_structure.py
new file mode 100644
index 000000000..de7f4865f
--- /dev/null
+++ b/sharpy/cases/templates/fuselage_wing_configuration/fwc_structure.py
@@ -0,0 +1,252 @@
+import numpy as np
+import h5py as h5
+
+
+class FWC_Structure:
+
+    """
+        FWC_Structure contains all attributes to define the beam geometriy 
+        and discretisation and makes it accessible for SHARPy.
+    """
+    def __init__(self, case_name, case_route, **kwargs):
+        self.n_elem_multiplier = kwargs.get('n_elem_multiplier', 1)
+
+        self.route = case_route
+        self.case_name = case_name
+
+        self.n_node_elem = 3
+        self.n_surfaces = 2
+        self.n_material = 2 # wing + fuselage
+        
+        self.half_wingspan = kwargs.get('half_wingspan', 2)
+        self.fuselage_length = kwargs.get('fuselage_length', 10)
+        self.offset_nose_wing_beam = kwargs.get('offset_nose_wing', self.fuselage_length/2)
+        self.vertical_wing_position = kwargs.get('vertical_wing_position', 0.)
+        self.n_elem_per_wing = kwargs.get('n_elem_per_wing', 10)
+        self.n_elem_fuselage = kwargs.get('n_elem_fuselage', 10)
+
+        self.sigma = kwargs.get('sigma', 10)
+        self.sigma_fuselage = kwargs.get('sigma_fuselage', 100.)
+
+
+
+        self.enforce_uniform_fuselage_discretisation = kwargs.get(
+            'enforce_uniform_fuselage_discretisation', False
+        )
+        self.fuselage_discretisation = kwargs.get('fuselage_discretisation', 'uniform')
+
+        self.thrust = 0.
+
+    def generate(self):
+        """
+            Function to set up all necessary parameter inputs to define the geometry and discretisation
+            of the beam. Finally, informations are saved to an .h5 file serving as an input
+            file for SHARPy.
+        """
+        self.set_element_and_nodes()
+        self.initialize_parameters()
+
+        self.set_stiffness_and_mass_propoerties()
+        self.set_beam_properties_right_wing()
+        self.mirror_wing_beam()                
+        self.app_forces[0] = [0, self.thrust, 0, 0, 0, 0]
+        
+        self.set_beam_properties_fuselage()
+        self.write_input_file()
+
+    def set_element_and_nodes(self):
+        """
+            Based on the specified number of elements of the wing and fuselage, the number of nodes of each 
+            component and the total number of elements and nodes are defined here.
+        """
+
+        self.n_node_fuselage = self.n_elem_fuselage*(self.n_node_elem - 1)
+        self.n_node_right_wing = self.n_elem_per_wing*(self.n_node_elem - 1) + 1
+        # the left wing beam has one node less than the right one, since they shares the center node
+        self.n_node_left_wing = self.n_node_right_wing - 1 
+
+        self.n_node_wing_total = self.n_node_right_wing + self.n_node_left_wing
+        self.n_node_fuselage_tail = self.n_node_wing_total + self.n_node_fuselage
+        self.n_node = self.n_node_fuselage + self.n_node_wing_total
+        self.n_elem =self.n_elem_fuselage + 2 * self.n_elem_per_wing
+
+        # vertical wing offset requires a connection element between wing and fuselage beam
+        if not self.vertical_wing_position == 0:
+            self.n_elem += 1
+            self.n_node += 1
+
+    def initialize_parameters(self):
+        self.x = np.zeros((self.n_node, ))
+        self.y = np.zeros((self.n_node, ))
+        self.z = np.zeros((self.n_node, ))
+
+        self.frame_of_reference_delta = np.zeros((self.n_elem, self.n_node_elem, 3))
+        self.conn = np.zeros((self.n_elem, self.n_node_elem), dtype=int)
+
+        self.beam_number = np.zeros((self.n_elem, ), dtype=int)
+        self.boundary_conditions = np.zeros((self.n_node, ), dtype=int)
+        self.structural_twist = np.zeros((self.n_elem, self.n_node_elem))
+
+        self.app_forces = np.zeros((self.n_node, 6))
+
+
+        self.stiffness = np.zeros((self.n_material, 6, 6))
+        self.mass = np.zeros((self.n_material, 6, 6))
+        self.elem_stiffness = np.zeros((self.n_elem, ), dtype=int)
+        self.mass = np.zeros((self.n_material, 6, 6))
+        self.elem_mass = np.zeros((self.n_elem, ), dtype=int)
+
+    def set_stiffness_and_mass_propoerties(self):
+        # Define aeroelastic properties
+        ea = 1e7
+        ga = 1e5
+        gj = 1e4
+        eiy = 2e4
+        eiz = 4e6
+        m_bar_main = 0.75
+        j_bar_main = 0.075
+
+        m_bar_fuselage = 0.3*1.5
+        j_bar_fuselage = 0.08
+
+        base_stiffness_main = self.sigma*np.diag([ea, ga, ga, gj, eiy, eiz])
+        base_stiffness_fuselage = base_stiffness_main.copy()*self.sigma_fuselage
+        base_stiffness_fuselage[4, 4] = base_stiffness_fuselage[5, 5]
+
+        self.stiffness[0, ...] = base_stiffness_main
+        self.stiffness[1, ...] = base_stiffness_fuselage
+
+        self.mass[0, ...] = self.generate_mass_matrix(m_bar_main, j_bar_main)
+        self.mass[1, ...] = self.generate_mass_matrix(m_bar_fuselage, j_bar_fuselage)
+        
+    def generate_mass_matrix(self, m_bar, j_bar):
+        return np.diag([m_bar, m_bar, m_bar, 
+                j_bar, 0.5*j_bar, 0.5*j_bar])
+
+    def set_beam_properties_right_wing(self):
+        """
+            Defines all necessary parameters to define the beam including node coordinate,
+            elements with their associated nodes, frame of reference delta, boundary conditions
+            such as reference node and free tips, stiffness and mass property ID for each element,
+            and twist.
+        """     
+        self.y[:self.n_node_right_wing] = np.linspace(0, self.half_wingspan, self.n_node_right_wing)
+        self.z[:self.n_node_right_wing] += self.vertical_wing_position
+        for ielem in range(self.n_elem_per_wing):
+            self.conn[ielem, :] = ((np.ones((3, )) * ielem * (self.n_node_elem - 1)) +
+                                [0, 2, 1])               
+            for ilocalnode in range(self.n_node_elem):
+                self.frame_of_reference_delta[ielem, ilocalnode, :] = [-1.0, 0.0, 0.0]  
+    
+        self.boundary_conditions[0] = 1
+        self.boundary_conditions[self.n_node_right_wing-1] = -1 # free tip
+
+    def mirror_wing_beam(self):
+        """
+            Mirrors the parameters from the beam representing the right free-flying wing
+            for the left one.
+        """
+
+        self.x[self.n_node_right_wing:self.n_node_wing_total]  = self.x[1:self.n_node_right_wing]
+        self.y[self.n_node_right_wing:self.n_node_wing_total]  = - self.y[1:self.n_node_right_wing]
+        self.z[self.n_node_right_wing:self.n_node_wing_total]  = self.z[1:self.n_node_right_wing]
+        self.frame_of_reference_delta[self.n_elem_per_wing:2*self.n_elem_per_wing, :, :] = self.frame_of_reference_delta[:self.n_elem_per_wing, :, :] * (-1)
+        self.elem_stiffness[self.n_elem_per_wing:2*self.n_elem_per_wing] = self.elem_stiffness[:self.n_elem_per_wing]
+        self.elem_mass[self.n_elem_per_wing:2*self.n_elem_per_wing] = self.elem_mass[:self.n_elem_per_wing]
+
+        self.beam_number[self.n_elem_per_wing:2*self.n_elem_per_wing] = 1
+        self.boundary_conditions[self.n_node_wing_total-1] = -1 # free tip
+        self.conn[self.n_elem_per_wing:2*self.n_elem_per_wing, :] = self.conn[:self.n_elem_per_wing, :] + self.n_node_right_wing - 1
+        self.conn[self.n_elem_per_wing, 0] = 0
+
+    def set_x_coordinate_fuselage(self):
+        if self.vertical_wing_position == 0:
+            n_nodes_fuselage = self.n_node_fuselage + 1
+        else:
+            n_nodes_fuselage = self.n_node_fuselage
+        if self.fuselage_discretisation == 'uniform':
+            x_coord_fuselage = np.linspace(0, self.fuselage_length, n_nodes_fuselage) - self.offset_nose_wing_beam
+        elif self.fuselage_discretisation == '1-cosine':
+            x_coord_fuselage = np.linspace(0, 1, n_nodes_fuselage) 
+            x_coord_fuselage =  0.5*(1.0 - np.cos(x_coord_fuselage*np.pi))
+            x_coord_fuselage *= self.fuselage_length
+            x_coord_fuselage -= self.offset_nose_wing_beam
+        else:
+            raise "ERROR Specified fuselage discretisation '{}' unknown".format(self.fuselage_discretisation)
+        self.idx_junction = self.find_index_of_closest_entry(x_coord_fuselage, self.x[0])
+        
+        self.idx_junction_global = self.idx_junction + self.n_node_wing_total
+        if self.vertical_wing_position == 0:
+            if self.enforce_uniform_fuselage_discretisation:
+                self.x[:self.n_node_wing_total] += x_coord_fuselage[self.idx_junction]
+            x_coord_fuselage = np.delete(x_coord_fuselage, self.idx_junction)
+        self.x[self.n_node_wing_total:self.n_node_fuselage_tail] = x_coord_fuselage
+
+    def adjust_fuselage_connectivities(self):
+        idx_in_conn = np.where(self.conn ==  self.idx_junction_global)
+        self.conn[idx_in_conn[0][0]+1:, :] -= 1
+
+        if idx_in_conn[0][0] == 2:
+            # if middle node, correct end node of element 
+            self.conn[idx_in_conn[0][0], 1] -= 1
+        for i_match in range(np.shape(idx_in_conn)[1]):
+            #several matches possible if junction node is not middle node
+            self.conn[idx_in_conn[0][i_match], idx_in_conn[1][i_match]] = 0
+
+    def add_additional_element_for_low_wing(self):
+        self.x[-1] = self.x[0]
+        self.y[-1] = self.y[0]
+        self.z[-1] = self.vertical_wing_position / 2
+        self.conn[-1, 0] = 0
+        self.conn[-1, 1] = self.idx_junction_global
+        self.conn[-1, 2] = self.n_node - 1
+        self.elem_stiffness[-1] = 1
+        self.elem_mass[-1] = 1
+        self.beam_number[-1] = 3
+        
+        
+    def set_beam_properties_fuselage(self):
+        self.set_x_coordinate_fuselage()
+        self.beam_number[self.n_elem_per_wing*2:] = 2
+        for ielem in range(self.n_elem_per_wing * 2,self.n_elem):
+            self.conn[ielem, :] = ((np.ones((3, )) * (ielem-self.n_elem_per_wing * 2) * (self.n_node_elem - 1)) +
+                                [0, 2, 1])  + self.n_node_wing_total
+            for ilocalnode in range(self.n_node_elem):
+                self.frame_of_reference_delta[ielem, ilocalnode, :] = [0.0, 1.0, 0.0]  
+        self.elem_stiffness[self.n_elem_per_wing*2:] = 1
+        self.elem_mass[self.n_elem_per_wing*2:] = 1
+        
+        if self.vertical_wing_position == 0:
+            self.adjust_fuselage_connectivities()
+        else:
+            self.add_additional_element_for_low_wing()
+        self.boundary_conditions[self.n_node_wing_total] = -1 # fuselage nose
+        self.boundary_conditions[self.n_node_wing_total + self.n_node_fuselage - 1] = -1 # fuselage tail
+
+    def find_index_of_closest_entry(self, array_values, target_value):
+        return np.argmin(np.abs(array_values - target_value))
+
+    def set_thrust(self, value):
+        self.thrust = value
+
+    def write_input_file(self):
+        """
+            Writes previously defined parameters to an .h5 file which serves later as an 
+            input file for SHARPy.
+        """                
+        with h5.File(self.route + '/' + self.case_name + '.fem.h5', 'a') as h5file:
+            h5file.create_dataset('coordinates', data=np.column_stack((self.x, self.y, self.z)))
+            h5file.create_dataset('connectivities', data=self.conn)
+            h5file.create_dataset('num_node_elem', data=self.n_node_elem)
+            h5file.create_dataset('num_node', data=self.n_node)
+            h5file.create_dataset('num_elem', data=self.n_elem)
+            h5file.create_dataset('stiffness_db', data=self.stiffness)
+            h5file.create_dataset('elem_stiffness', data=self.elem_stiffness)
+            h5file.create_dataset('mass_db', data=self.mass)
+            h5file.create_dataset('elem_mass', data=self.elem_mass)
+            h5file.create_dataset('frame_of_reference_delta', data=self.frame_of_reference_delta)
+            h5file.create_dataset('boundary_conditions', data=self.boundary_conditions)
+            h5file.create_dataset('beam_number', data=self.beam_number)
+
+            h5file.create_dataset('structural_twist', data=self.structural_twist)
+            h5file.create_dataset('app_forces', data=self.app_forces)
diff --git a/sharpy/generators/polaraeroforces.py b/sharpy/generators/polaraeroforces.py
index b17e59ff0..b94c368d7 100644
--- a/sharpy/generators/polaraeroforces.py
+++ b/sharpy/generators/polaraeroforces.py
@@ -155,7 +155,7 @@ def generate(self, **params):
         
         list_aoa_induced = []
 
-        aero_dict = aerogrid.aero_dict
+        data_dict = aerogrid.data_dict
         if aerogrid.polars is None:
             return struct_forces
         new_struct_forces = np.zeros_like(struct_forces)
@@ -168,9 +168,9 @@ def generate(self, **params):
 
         for inode in range(nnode):
             new_struct_forces[inode, :] = struct_forces[inode, :].copy()
-            if aero_dict['aero_node'][inode]:
+            if data_dict['aero_node'][inode]:
                 ielem, inode_in_elem = structure.node_master_elem[inode]
-                iairfoil = aero_dict['airfoil_distribution'][ielem, inode_in_elem]
+                iairfoil = data_dict['airfoil_distribution'][ielem, inode_in_elem]
                 isurf = aerogrid.struct2aero_mapping[inode][0]['i_surf']
                 if isurf not in self.settings['skip_surfaces']:
                     i_n = aerogrid.struct2aero_mapping[inode][0]['i_n']
@@ -179,7 +179,7 @@ def generate(self, **params):
                     cgb = np.dot(cga, cab)
 
                     if not self.cd_from_cl:
-                        airfoil = str(aero_dict['airfoil_distribution'][ielem, inode_in_elem])
+                        airfoil = str(data_dict['airfoil_distribution'][ielem, inode_in_elem])
                         aoa_0cl = self.list_aoa_cl0[int(airfoil)]
 
                     # computing surface area of panels contributing to force
@@ -317,10 +317,10 @@ def check_for_special_cases(self, aerogrid):
         # check if outboard node of aerosurface
         self.flag_shared_node_by_surfaces = np.zeros((self.n_node,1))
         for inode in range(self.n_node):
-            if aerogrid.aero_dict['aero_node'][inode]:
+            if aerogrid.data_dict['aero_node'][inode]:
                 i_n = aerogrid.struct2aero_mapping[inode][0]['i_n']                
                 isurf = aerogrid.struct2aero_mapping[inode][0]['i_surf']
-                N = aerogrid.aero_dimensions[isurf, 1]
+                N = aerogrid.dimensions[isurf, 1]
                 if i_n in [0, N]:
                     if len(aerogrid.struct2aero_mapping[inode]) > 1:
                         self.flag_shared_node_by_surfaces[inode] = 1
@@ -350,9 +350,9 @@ def compute_aoa_cl0_from_airfoil_data(self, aerogrid):
         """
          Computes the angle of attack for which zero lift is achieved for every airfoil
         """
-        self.list_aoa_cl0 = np.zeros((len(aerogrid.aero_dict['airfoils']),1))
-        for i, airfoil in enumerate(aerogrid.aero_dict['airfoils']):
-            airfoil_coords = aerogrid.aero_dict['airfoils'][airfoil]
+        self.list_aoa_cl0 = np.zeros((len(aerogrid.data_dict['airfoils']),1))
+        for i, airfoil in enumerate(aerogrid.data_dict['airfoils']):
+            airfoil_coords = aerogrid.data_dict['airfoils'][airfoil]
             self.list_aoa_cl0[i] = get_aoacl0_from_camber(airfoil_coords[:, 0], airfoil_coords[:, 1])
 
 @generator_interface.generator
@@ -401,11 +401,11 @@ def generate(self, **params):
 
         n_node = self.structure.num_node
         n_elem = self.structure.num_elem
-        aero_dict = self.aero.aero_dict
+        data_dict = self.aero.data_dict
         new_struct_forces = np.zeros_like(struct_forces)
 
         # load airfoil efficiency (if it exists); else set to one (to avoid multiple ifs in the loops)
-        airfoil_efficiency = aero_dict['airfoil_efficiency']
+        airfoil_efficiency = data_dict['airfoil_efficiency']
         # force efficiency dimensions [n_elem, n_node_elem, 2, [fx, fy, fz]] - all defined in B frame
         force_efficiency = np.zeros((n_elem, 3, 2, 3))
         force_efficiency[:, :, 0, :] = 1.
diff --git a/sharpy/linear/assembler/lincontrolsurfacedeflector.py b/sharpy/linear/assembler/lincontrolsurfacedeflector.py
index 11dbb84fd..2f99fda4b 100644
--- a/sharpy/linear/assembler/lincontrolsurfacedeflector.py
+++ b/sharpy/linear/assembler/lincontrolsurfacedeflector.py
@@ -82,7 +82,7 @@ def generate(self):
         tsstruct0 = self.tsstruct0
 
         # Find the vertices corresponding to a control surface from beam coordinates to aerogrid
-        aero_dict = aero.aero_dict
+        data_dict = aero.data_dict
         n_surf = tsaero0.n_surf
         n_control_surfaces = self.n_control_surfaces
 
@@ -120,29 +120,29 @@ def generate(self):
 
                         # Although a node may be part of 2 aerodynamic surfaces, we need to ensure that the current
                         # element for the given node is indeed part of that surface.
-                        elems_in_surf = np.where(aero_dict['surface_distribution'] == i_surf)[0]
+                        elems_in_surf = np.where(data_dict['surface_distribution'] == i_surf)[0]
                         if i_elem not in elems_in_surf:
                             continue
 
                         # Surface panelling
-                        M = aero.aero_dimensions[i_surf][0]
-                        N = aero.aero_dimensions[i_surf][1]
+                        M = aero.dimensions[i_surf][0]
+                        N = aero.dimensions[i_surf][1]
 
                         K_zeta_start = 3 * sum(linuvlm.MS.KKzeta[:i_surf])
                         shape_zeta = (3, M + 1, N + 1)
 
-                        i_control_surface = aero_dict['control_surface'][i_elem, i_local_node]
+                        i_control_surface = data_dict['control_surface'][i_elem, i_local_node]
                         if i_control_surface >= 0:
                             if not with_control_surface:
                                 i_start_of_cs = i_node_span.copy()
                                 with_control_surface = True
 
-                            control_surface_chord = aero_dict['control_surface_chord'][i_control_surface]
+                            control_surface_chord = data_dict['control_surface_chord'][i_control_surface]
 
                             try:
                                 control_surface_hinge_coord = \
-                                    aero_dict['control_surface_hinge_coord'][i_control_surface] * \
-                                    aero_dict['chord'][i_elem, i_local_node]
+                                    data_dict['control_surface_hinge_coord'][i_control_surface] * \
+                                    data_dict['chord'][i_elem, i_local_node]
                             except KeyError:
                                 control_surface_hinge_coord = None
 
diff --git a/sharpy/linear/assembler/linearaeroelastic.py b/sharpy/linear/assembler/linearaeroelastic.py
index 752654cd1..0dfdfd1ae 100644
--- a/sharpy/linear/assembler/linearaeroelastic.py
+++ b/sharpy/linear/assembler/linearaeroelastic.py
@@ -680,8 +680,8 @@ def get_gebm2uvlm_gains(self, data):
                 # print('%.2d,%.2d'%(nn,ss))
 
                 # surface panelling
-                M = aero.aero_dimensions[ss][0]
-                N = aero.aero_dimensions[ss][1]
+                M = aero.dimensions[ss][0]
+                N = aero.dimensions[ss][1]
 
                 Kzeta_start = 3 * sum(self.uvlm.sys.MS.KKzeta[:ss])
                 shape_zeta = (3, M + 1, N + 1)
diff --git a/sharpy/linear/assembler/lineargustassembler.py b/sharpy/linear/assembler/lineargustassembler.py
index 914ceceeb..72b6878ca 100644
--- a/sharpy/linear/assembler/lineargustassembler.py
+++ b/sharpy/linear/assembler/lineargustassembler.py
@@ -250,7 +250,7 @@ def assemble(self):
         c_i = np.zeros((3 * Kzeta, N))
         for i_surf in range(self.aero.n_surf):
 
-            M_surf, N_surf = self.aero.aero_dimensions[i_surf]
+            M_surf, N_surf = self.aero.dimensions[i_surf]
             Kzeta_start = 3 * sum(self.KKzeta[:i_surf])  # number of coordinates up to current surface
             shape_zeta = (3, M_surf + 1, N_surf + 1)
 
@@ -341,7 +341,7 @@ def assemble(self):
         gust_d = np.zeros((gust_c.shape[0], gust_b.shape[1]))
         for i_surf in range(self.aero.n_surf):
 
-            M_surf, N_surf = self.aero.aero_dimensions[i_surf]
+            M_surf, N_surf = self.aero.dimensions[i_surf]
             Kzeta_start = 3 * sum(self.KKzeta[:i_surf])  # number of coordinates up to current surface
             shape_zeta = (3, M_surf + 1, N_surf + 1)
 
diff --git a/sharpy/linear/assembler/linearuvlm.py b/sharpy/linear/assembler/linearuvlm.py
index b9b20e277..59df15366 100644
--- a/sharpy/linear/assembler/linearuvlm.py
+++ b/sharpy/linear/assembler/linearuvlm.py
@@ -425,8 +425,8 @@ def unpack_ss_vector(self, data, x_n, u_aero, aero_tstep, track_body=False, stat
         for i_surf in range(aero_tstep.n_surf):
             # Tuple with dimensions of the aerogrid zeta, which is the same shape for forces
             dimensions = aero_tstep.zeta[i_surf].shape
-            dimensions_gamma = data.aero.aero_dimensions[i_surf]
-            dimensions_wake = data.aero.aero_dimensions_star[i_surf]
+            dimensions_gamma = data.aero.dimensions[i_surf]
+            dimensions_wake = data.aero.dimensions_star[i_surf]
 
             # Number of entries in zeta
             points_in_surface = aero_tstep.zeta[i_surf].size
diff --git a/sharpy/postproc/aeroforcescalculator.py b/sharpy/postproc/aeroforcescalculator.py
index fce5468b8..219b0f783 100644
--- a/sharpy/postproc/aeroforcescalculator.py
+++ b/sharpy/postproc/aeroforcescalculator.py
@@ -6,6 +6,7 @@
 import sharpy.utils.settings as settings_utils
 import sharpy.utils.algebra as algebra
 import sharpy.aero.utils.mapping as mapping
+import warnings
 
 
 @solver
@@ -38,6 +39,14 @@ class AeroForcesCalculator(BaseSolver):
     settings_types['coefficients'] = 'bool'
     settings_description['coefficients'] = 'Calculate aerodynamic coefficients'
 
+    settings_default['lifting_surfaces'] = True
+    settings_types['lifting_surfaces'] = 'bool'
+    settings_description['lifting_surfaces'] = 'Includes aerodynamic forces from lifting surfaces'
+
+    settings_default['nonlifting_body'] = False
+    settings_types['nonlifting_body'] = 'bool'
+    settings_description['nonlifting_body'] = 'Includes aerodynamic forces from nonlifting bodies'
+
     settings_types['q_ref'] = 'float'
     settings_default['q_ref'] = 1
     settings_description['q_ref'] = 'Reference dynamic pressure'
@@ -67,6 +76,8 @@ def __init__(self):
         self.caller = None
 
         self.table = None
+        self.rot = None
+        self.moment_reference_location = np.array([0., 0., 0.])
 
     def initialise(self, data, custom_settings=None, caller=None, restart=False):
         self.data = data
@@ -107,38 +118,43 @@ def run(self, **kwargs):
         if self.settings['write_text_file']:
             self.file_output(self.settings['text_file_name'])
         return self.data
-
+ 
     def calculate_forces(self, ts):
-        rot = algebra.quat2rotation(self.data.structure.timestep_info[ts].quat)
-
         # Forces per surface in G frame
+        self.rot = algebra.quat2rotation(self.data.structure.timestep_info[ts].quat)
         force = self.data.aero.timestep_info[ts].forces
         unsteady_force = self.data.aero.timestep_info[ts].dynamic_forces
-        n_surf = len(force)
-        for i_surf in range(n_surf):
-            total_steady_force = np.zeros((3,))
-            total_unsteady_force = np.zeros((3,))
-            _, n_rows, n_cols = force[i_surf].shape
-            for i_m in range(n_rows):
-                for i_n in range(n_cols):
-                    total_steady_force += force[i_surf][0:3, i_m, i_n]
-                    total_unsteady_force += unsteady_force[i_surf][0:3, i_m, i_n]
-            self.data.aero.timestep_info[ts].inertial_steady_forces[i_surf, 0:3] = total_steady_force
-            self.data.aero.timestep_info[ts].inertial_unsteady_forces[i_surf, 0:3] = total_unsteady_force
-            self.data.aero.timestep_info[ts].body_steady_forces[i_surf, 0:3] = np.dot(rot.T, total_steady_force)
-            self.data.aero.timestep_info[ts].body_unsteady_forces[i_surf, 0:3] = np.dot(rot.T, total_unsteady_force)
-
-        # Forces expressed in the beam degrees of freedom
+        for i_surf in range(self.data.aero.n_surf):
+            (
+            self.data.aero.timestep_info[ts].inertial_steady_forces[i_surf, 0:3], 
+            self.data.aero.timestep_info[ts].inertial_unsteady_forces[i_surf, 0:3],
+            self.data.aero.timestep_info[ts].body_steady_forces[i_surf, 0:3],
+            self.data.aero.timestep_info[ts].body_unsteady_forces[i_surf, 0:3]
+            ) = self.calculate_forces_for_isurf_in_g_frame(force[i_surf], unsteady_force=unsteady_force[i_surf])
+        
+        print(self.settings["nonlifting_body"])
+        if self.settings["nonlifting_body"]:
+            print(self.data.nonlifting_body.n_surf)
+            for i_surf in range(self.data.nonlifting_body.n_surf):
+                print(i_surf)
+                (
+                self.data.nonlifting_body.timestep_info[ts].inertial_steady_forces[i_surf, 0:3], 
+                self.data.nonlifting_body.timestep_info[ts].body_steady_forces[i_surf, 0:3],
+                ) = self.calculate_forces_for_isurf_in_g_frame(self.data.nonlifting_body.timestep_info[ts].forces[i_surf], nonlifting=True)
+        
+        # Convert to forces in B frame
         try:
             steady_forces_b = self.data.structure.timestep_info[ts].postproc_node['aero_steady_forces']
         except KeyError:
-            steady_forces_b = self.map_forces_beam_dof(ts, force)
+            if self.settings["nonlifting_body"]:
+                warnings.warn('Nonlifting forces are not considered in aero forces calculation since forces cannot not be retrieved from postproc node.')
+            steady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, force)
 
         try:
             unsteady_forces_b = self.data.structure.timestep_info[ts].postproc_node['aero_unsteady_forces']
         except KeyError:
-            unsteady_forces_b = self.map_forces_beam_dof(ts, unsteady_force)
-
+            unsteady_forces_b = self.map_forces_beam_dof(self.data.aero, ts, unsteady_force)
+        # Convert to forces in A frame
         steady_forces_a = self.data.structure.timestep_info[ts].nodal_b_for_2_a_for(steady_forces_b,
                                                                                     self.data.structure)
 
@@ -146,33 +162,50 @@ def calculate_forces(self, ts):
                                                                                       self.data.structure)
 
         # Express total forces in A frame
-        moment_reference_location = np.array([0., 0., 0.])
         self.data.aero.timestep_info[ts].total_steady_body_forces = \
             mapping.total_forces_moments(steady_forces_a,
                                          self.data.structure.timestep_info[ts].pos,
-                                         ref_pos=moment_reference_location)
+                                         ref_pos=self.moment_reference_location)
         self.data.aero.timestep_info[ts].total_unsteady_body_forces = \
             mapping.total_forces_moments(unsteady_forces_a,
                                          self.data.structure.timestep_info[ts].pos,
-                                         ref_pos=moment_reference_location)
+                                         ref_pos=self.moment_reference_location)
 
         # Express total forces in G frame
         self.data.aero.timestep_info[ts].total_steady_inertial_forces = \
-            np.block([[rot, np.zeros((3, 3))],
-                      [np.zeros((3, 3)), rot]]).dot(
+            np.block([[self.rot, np.zeros((3, 3))],
+                      [np.zeros((3, 3)), self.rot]]).dot(
                 self.data.aero.timestep_info[ts].total_steady_body_forces)
 
         self.data.aero.timestep_info[ts].total_unsteady_inertial_forces = \
-            np.block([[rot, np.zeros((3, 3))],
-                      [np.zeros((3, 3)), rot]]).dot(
+            np.block([[self.rot, np.zeros((3, 3))],
+                      [np.zeros((3, 3)), self.rot]]).dot(
                 self.data.aero.timestep_info[ts].total_unsteady_body_forces)
 
-    def map_forces_beam_dof(self, ts, force):
-        aero_tstep = self.data.aero.timestep_info[ts]
+    def calculate_forces_for_isurf_in_g_frame(self, force, unsteady_force = None, nonlifting = False):
+        """
+            Forces for a surface in G frame
+        """
+        # Forces per surface in G frame
+        total_steady_force = np.zeros((3,))
+        total_unsteady_force = np.zeros((3,))
+        _, n_rows, n_cols = force.shape
+        for i_m in range(n_rows):
+            for i_n in range(n_cols):
+                total_steady_force += force[0:3, i_m, i_n]
+                if not nonlifting:
+                    total_unsteady_force += unsteady_force[0:3, i_m, i_n]
+        if not nonlifting:
+            return total_steady_force, total_unsteady_force, np.dot(self.rot.T, total_steady_force), np.dot(self.rot.T, total_unsteady_force)
+        else:
+            return total_steady_force, np.dot(self.rot.T, total_steady_force)
+
+
+    def map_forces_beam_dof(self, aero_data, ts, force):
         struct_tstep = self.data.structure.timestep_info[ts]
         aero_forces_beam_dof = mapping.aero2struct_force_mapping(force,
-                                                                 self.data.aero.struct2aero_mapping,
-                                                                 aero_tstep.zeta,
+                                                                 aero_data.struct2aero_mapping,
+                                                                 aero_data.timestep_info[ts].zeta,
                                                                  struct_tstep.pos,
                                                                  struct_tstep.psi,
                                                                  None,
@@ -187,15 +220,18 @@ def calculate_coefficients(self, fx, fy, fz, mx, my, mz):
 
     def screen_output(self, ts):
         # print time step total aero forces
+        forces  = np.zeros(3)
+        moments  = np.zeros(3)
+
         aero_tstep = self.data.aero.timestep_info[ts]
-        fx, fy, fz = aero_tstep.total_steady_inertial_forces[:3] + aero_tstep.total_unsteady_inertial_forces[:3]
-        mx, my, mz = aero_tstep.total_steady_inertial_forces[3:] + aero_tstep.total_unsteady_inertial_forces[3:]
+        forces += aero_tstep.total_steady_inertial_forces[:3] + aero_tstep.total_unsteady_inertial_forces[:3]
+        moments += aero_tstep.total_steady_inertial_forces[3:] + aero_tstep.total_unsteady_inertial_forces[3:]
 
-        if self.settings['coefficients']:
-            Cfx, Cfy, Cfz, Cmx, Cmy, Cmz = self.calculate_coefficients(fx, fy, fz, mx, my, mz)
+        if self.settings['coefficients']: # TODO: Check if coefficients have to be computed differently for fuselages
+            Cfx, Cfy, Cfz, Cmx, Cmy, Cmz = self.calculate_coefficients(*forces, *moments)
             self.table.print_line([ts, Cfx, Cfy, Cfz, Cmx, Cmy, Cmz])
         else:
-            self.table.print_line([ts, fx, fy, fz, mx, my, mz])
+            self.table.print_line([ts, *forces, *moments])
 
     def file_output(self, filename):
         # assemble forces/moments matrix
diff --git a/sharpy/postproc/aerogridplot.py b/sharpy/postproc/aerogridplot.py
index ea925dfe7..f860459f3 100644
--- a/sharpy/postproc/aerogridplot.py
+++ b/sharpy/postproc/aerogridplot.py
@@ -65,6 +65,14 @@ class AerogridPlot(BaseSolver):
     settings_default['include_incidence_angle'] = False
     settings_description['include_incidence_angle'] = 'Include panel incidence angle'
 
+    settings_types['plot_nonlifting_surfaces'] = 'bool'
+    settings_default['plot_nonlifting_surfaces'] = False
+    settings_description['plot_nonlifting_surfaces'] = 'Plot nonlifting surfaces'
+    
+    settings_types['plot_lifting_surfaces'] = 'bool'
+    settings_default['plot_lifting_surfaces'] = False
+    settings_description['plot_lifting_surfaces'] = 'Plot nonlifting surfaces'
+
     settings_types['num_cores'] = 'int'
     settings_default['num_cores'] = 1
     settings_description['num_cores'] = 'Number of cores used to compute velocities/angles'
@@ -91,6 +99,7 @@ def __init__(self):
         self.folder = None
         self.body_filename = ''
         self.wake_filename = ''
+        self.nonlifting_filename = ''
         self.ts_max = 0
         self.caller = None
 
@@ -114,6 +123,10 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
                               self.settings['name_prefix'] +
                               'wake_' +
                               self.data.settings['SHARPy']['case'])
+        self.nonlifting_filename = (self.folder +
+                                    self.settings['name_prefix'] +
+                                    'nonlifting_' +
+                                    self.data.settings['SHARPy']['case'])
         self.caller = caller
 
     def run(self, **kwargs):
@@ -125,16 +138,18 @@ def run(self, **kwargs):
             for self.ts in range(self.ts_max):
                 if self.data.structure.timestep_info[self.ts] is not None:
                     self.plot_body()
-                    if self.settings['save_wake']:
-                        self.plot_wake()
+                    self.plot_wake()
+                    if self.settings['plot_nonlifting_surfaces']:
+                        self.plot_nonlifting_surfaces()
             cout.cout_wrap('...Finished', 1)
         elif (self.data.ts % self.settings['stride'] == 0):
             aero_tsteps = len(self.data.aero.timestep_info) - 1
             struct_tsteps = len(self.data.structure.timestep_info) - 1
             self.ts = np.max((aero_tsteps, struct_tsteps))
             self.plot_body()
-            if self.settings['save_wake']:
-                self.plot_wake()
+            self.plot_wake()
+            if self.settings['plot_nonlifting_surfaces']:
+                self.plot_nonlifting_surfaces()
         return self.data
 
     def plot_body(self):
@@ -333,3 +348,105 @@ def plot_wake(self):
             ug.point_data.scalars = np.arange(0, coords.shape[0])
             ug.point_data.scalars.name = 'n_id'
             write_data(ug, filename)
+
+    def plot_nonlifting_surfaces(self):
+        nonlifting_tstep = self.data.nonlifting_body.timestep_info[self.ts]
+        struct_tstep = self.data.structure.timestep_info[self.ts]
+
+        for i_surf in range(nonlifting_tstep.n_surf):
+            filename = (self.nonlifting_filename +
+                        '_' +
+                        '%02u_' % i_surf +
+                        '%06u' % self.ts+
+                        '.vtu')
+
+            dims = nonlifting_tstep.dimensions[i_surf, :]
+            point_data_dim = (dims[0]+1)*(dims[1]+1)  # + (dims_star[0]+1)*(dims_star[1]+1)
+            panel_data_dim = (dims[0])*(dims[1])  # + (dims_star[0])*(dims_star[1])
+
+            coords = np.zeros((point_data_dim, 3))
+            conn = []
+            panel_id = np.zeros((panel_data_dim,), dtype=int)
+            panel_surf_id = np.zeros((panel_data_dim,), dtype=int)
+            panel_sigma = np.zeros((panel_data_dim,))
+            normal = np.zeros((panel_data_dim, 3))
+            point_struct_id = np.zeros((point_data_dim,), dtype=int)
+            point_cf = np.zeros((point_data_dim, 3))
+            u_inf = np.zeros((point_data_dim, 3))
+            counter = -1
+
+            # coordinates of corners
+            for i_n in range(dims[1]+1):
+                for i_m in range(dims[0]+1):
+                    counter += 1
+                    coords[counter, :] = nonlifting_tstep.zeta[i_surf][:, i_m, i_n]
+                    # TODO: include those for nonlifting body (are they different for nonlifting coordinates?)
+                    if self.settings['include_rbm']:
+                        coords[counter, :] += struct_tstep.for_pos[0:3]
+                    if self.settings['include_forward_motion']:
+                        coords[counter, 0] -= self.settings['dt']*self.ts*self.settings['u_inf']
+            counter = -1
+            node_counter = -1
+            for i_n in range(dims[1] + 1):
+                global_counter = self.data.nonlifting_body.aero2struct_mapping[i_surf][i_n]
+                for i_m in range(dims[0] + 1):
+                    node_counter += 1
+                    # point data
+                    point_struct_id[node_counter] = global_counter
+                    point_cf[node_counter, :] = nonlifting_tstep.forces[i_surf][0:3, i_m, i_n]
+                    try:
+                        u_inf[node_counter, :] = nonlifting_tstep.u_ext[i_surf][0:3, i_m, i_n]
+                    except AttributeError:
+                        pass
+                    if i_n < dims[1] and i_m < dims[0]:
+                        counter += 1
+                    else:
+                        continue
+
+                    conn.append([node_counter + 0,
+                                 node_counter + 1,
+                                 node_counter + dims[0]+2,
+                                 node_counter + dims[0]+1])
+                    # cell data
+                    normal[counter, :] = nonlifting_tstep.normals[i_surf][:, i_m, i_n]
+                    panel_id[counter] = counter
+                    panel_surf_id[counter] = i_surf
+                    panel_sigma[counter] = nonlifting_tstep.sigma[i_surf][i_m, i_n]
+
+            ug = tvtk.UnstructuredGrid(points=coords)
+            ug.set_cells(tvtk.Quad().cell_type, conn)
+            ug.cell_data.scalars = panel_id
+            ug.cell_data.scalars.name = 'panel_n_id'
+            ug.cell_data.add_array(panel_surf_id)
+            ug.cell_data.get_array(1).name = 'panel_surface_id'
+            ug.cell_data.add_array(panel_sigma)
+            ug.cell_data.get_array(2).name = 'panel_sigma'
+            ug.cell_data.vectors = normal
+            ug.cell_data.vectors.name = 'panel_normal'
+            ug.point_data.scalars = np.arange(0, coords.shape[0])
+            ug.point_data.scalars.name = 'n_id'
+            ug.point_data.add_array(point_struct_id)
+            ug.point_data.get_array(1).name = 'point_struct_id'
+            ug.point_data.add_array(point_cf)
+            ug.point_data.get_array(2).name = 'point_steady_force'
+            ug.point_data.add_array(u_inf)
+            ug.point_data.get_array(3).name = 'u_inf'
+            write_data(ug, filename)
+
+        def write_paraview_data(self, coords, conn, panel_id, list_cell_parameters, list_cell_names, list_point_parameters, list_point_names, filename):
+            ug = tvtk.UnstructuredGrid(points=coords)
+            
+            ug.set_cells(tvtk.Quad().cell_type, conn)
+            ug.cell_data.scalars = panel_id
+            ug.cell_data.scalars.name = 'panel_n_id'
+            for counter in range(len(list_cell_parameters)):
+                ug.cell_data.add_array(list_cell_parameters[counter])
+                ug.cell_data.get_array(counter+1).name = list_cell_names[counter]
+                
+            ug.point_data.scalars = np.arange(0, coords.shape[0])
+            ug.point_data.scalars.name = 'n_id'
+            for counter in range(len(list_point_parameters)):
+                ug.point_data.add_array(list_point_parameters[counter])
+                ug.point_data.get_array(counter+1).name = list_point_names[counter]
+            
+            write_data(ug, filename)
diff --git a/sharpy/postproc/liftdistribution.py b/sharpy/postproc/liftdistribution.py
index d3e2c7cd0..d0cf7e79d 100644
--- a/sharpy/postproc/liftdistribution.py
+++ b/sharpy/postproc/liftdistribution.py
@@ -2,21 +2,19 @@
 
 import numpy as np
 
-import sharpy.utils.cout_utils as cout
-import sharpy.utils.algebra as algebra
 from sharpy.utils.solver_interface import solver, BaseSolver
 import sharpy.utils.settings as settings_utils
-from sharpy.utils.datastructures import init_matrix_structure, standalone_ctypes_pointer
 import sharpy.aero.utils.mapping as mapping
+import sharpy.utils.algebra as algebra
 import sharpy.aero.utils.utils as aeroutils
 
 
 @solver
 class LiftDistribution(BaseSolver):
     """LiftDistribution
-
+    
     Calculates and exports the lift distribution on lifting surfaces
-
+    
     """
     solver_id = 'LiftDistribution'
     solver_classification = 'post-processor'
@@ -26,7 +24,7 @@ class LiftDistribution(BaseSolver):
     settings_description = dict()
 
     settings_types['text_file_name'] = 'str'
-    settings_default['text_file_name'] = 'lift_distribution.csv'
+    settings_default['text_file_name'] = 'liftdistribution'
     settings_description['text_file_name'] = 'Text file name'
 
     settings_default['coefficients'] = True
@@ -46,30 +44,18 @@ def __init__(self):
         self.folder = None
         self.caller = None
 
-    def initialise(self, data, custom_settings=None, caller=None, restart=False):
+    def initialise(self, data, custom_settings=None, restart=False, caller=None):
         self.data = data
-        if custom_settings is None:
-            self.settings = data.settings[self.solver_id]
-        else:
-            self.settings = custom_settings
+        self.settings = data.settings[self.solver_id]
         settings_utils.to_custom_types(self.settings, self.settings_types, self.settings_default)
-        self.ts_max = len(self.data.structure.timestep_info)
         self.caller = caller
-        self.folder = data.output_folder
+        self.folder = data.output_folder + '/liftdistribution/'
         if not os.path.exists(self.folder):
             os.makedirs(self.folder)
 
     def run(self, **kwargs):
-
-        online = settings_utils.set_value_or_default(kwargs, 'online', False)
-
-        if not online:
-            for self.ts in range(self.ts_max):
-                self.lift_distribution()
-            cout.cout_wrap('...Finished', 1)
-        else:
-            self.ts = len(self.data.structure.timestep_info) - 1
-            self.lift_distribution()
+        self.lift_distribution(self.data.structure.timestep_info[self.data.ts],
+                               self.data.aero.timestep_info[self.data.ts])
         return self.data
 
     def lift_distribution(self, struct_tstep, aero_tstep):
@@ -83,27 +69,26 @@ def lift_distribution(self, struct_tstep, aero_tstep):
             self.data.structure.node_master_elem,
             self.data.structure.connectivities,
             struct_tstep.cag(),
-            self.data.aero.aero_dict)
+            self.data.aero.data_dict)
         # Prepare output matrix and file
         N_nodes = self.data.structure.num_node
-        numb_col = 4
-        header = "x,y,z,fz"
+        numb_col = 6
+        header = "x,y,z,fx,fy,fz"
         # get aero forces
-        lift_distribution = np.zeros((N_nodes, numb_col))
         # get rotation matrix
         cga = algebra.quat2rotation(struct_tstep.quat)
         if self.settings["coefficients"]:
             # TODO: add nondimensional spanwise column y/s
-            header += ", y/s, cl"
-            numb_col += 2
-            lift_distribution = np.concatenate((lift_distribution, np.zeros((N_nodes, 2))), axis=1)
+            header += ", cfx, cfy, cfz"
+            numb_col += 3
+        lift_distribution = np.zeros((N_nodes, numb_col))
 
         for inode in range(N_nodes):
-            if self.data.aero.aero_dict['aero_node'][inode]:
+            if self.data.aero.data_dict['aero_node'][inode]:
                 local_node = self.data.aero.struct2aero_mapping[inode][0]["i_n"]
                 ielem, inode_in_elem = self.data.structure.node_master_elem[inode]
                 i_surf = int(self.data.aero.surface_distribution[ielem])
-                # get c_gb
+                # get c_gb                
                 cab = algebra.crv2rotation(struct_tstep.psi[ielem, inode_in_elem, :])
                 cgb = np.dot(cga, cab)
                 # Get c_bs
@@ -116,23 +101,22 @@ def lift_distribution(self, struct_tstep, aero_tstep):
                 dir_span, span, dir_chord, chord = aeroutils.span_chord(local_node, aero_tstep.zeta[i_surf])
                 # Stability axes - projects forces in B onto S
                 c_bs = aeroutils.local_stability_axes(cgb.T.dot(dir_urel), cgb.T.dot(dir_chord))
-                lift_force = c_bs.T.dot(forces[inode, :3])[2]
+                aero_forces = c_bs.T.dot(forces[inode, :3])
                 # Store data in export matrix
-                lift_distribution[inode, 3] = lift_force
+                lift_distribution[inode, 3:6] = aero_forces
                 lift_distribution[inode, 2] = struct_tstep.pos[inode, 2]  # z
                 lift_distribution[inode, 1] = struct_tstep.pos[inode, 1]  # y
                 lift_distribution[inode, 0] = struct_tstep.pos[inode, 0]  # x
                 if self.settings["coefficients"]:
-                    # Get non-dimensional spanwise coordinate y/s
-                    lift_distribution[inode, 4] = lift_distribution[inode, 1]/span
                     # Get lift coefficient
-                    lift_distribution[inode, 5] = np.sign(lift_force) * np.linalg.norm(lift_force) \
-                                                  / (0.5 * self.settings['rho'] \
-                                                     * np.linalg.norm(urel) ** 2 * span * chord)  # strip_area[i_surf][local_node])
-                    # Check if shared nodes from different surfaces exist (e.g. two wings joining at symmetry plane)
-                    # Leads to error since panel area just donates for half the panel size while lift forces is summed up
-                    lift_distribution[inode, 5] /= len(self.data.aero.struct2aero_mapping[inode])
+                    for idim in range(3):
+                        lift_distribution[inode, 6+idim] = np.sign(aero_forces[idim]) * np.linalg.norm(aero_forces[idim]) \
+                                                    / (0.5 * self.settings['rho'] \
+                                                        * np.linalg.norm(urel) ** 2 * span * chord)  
+                        # Check if shared nodes from different surfaces exist (e.g. two wings joining at symmetry plane)
+                        # Leads to error since panel area just donates for half the panel size while lift forces is summed up
+                        lift_distribution[inode, 6+idim] /= len(self.data.aero.struct2aero_mapping[inode])
 
         # Export lift distribution data
-        np.savetxt(os.path.join(self.folder, self.settings['text_file_name']), lift_distribution,
+        np.savetxt(os.path.join(self.folder,  self.settings['text_file_name'] + '_ts{}'.format(str(self.data.ts)) + '.txt'), lift_distribution,
                    fmt='%10e,' * (numb_col - 1) + '%10e', delimiter=", ", header=header)
diff --git a/sharpy/postproc/savedata.py b/sharpy/postproc/savedata.py
index 083dbd4e8..87ccd2df5 100644
--- a/sharpy/postproc/savedata.py
+++ b/sharpy/postproc/savedata.py
@@ -40,6 +40,10 @@ class SaveData(BaseSolver):
     settings_default['save_aero'] = True
     settings_description['save_aero'] = 'Save aerodynamic classes.'
 
+    settings_types['save_nonlifting'] = 'bool'
+    settings_default['save_nonlifting'] = False
+    settings_description['save_nonlifting'] = 'Save aerodynamic classes.'
+
     settings_types['save_struct'] = 'bool'
     settings_default['save_struct'] = True
     settings_description['save_struct'] = 'Save structural classes.'
@@ -106,7 +110,6 @@ def __init__(self):
         self.folder = ''
         self.filename = ''
         self.filename_linear = ''
-        self.ts_max = 0
         self.caller = None
 
         ### specify which classes are saved as hdf5 group
@@ -143,7 +146,6 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
                                             'ct_zeta_star_list',
                                             'dynamic_input'])
 
-        self.ts_max = self.data.ts + 1
 
         # create folder for containing files if necessary
         self.folder = data.output_folder + '/savedata/'
@@ -180,6 +182,11 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
                     self.settings['skip_attr'].append('dist_to_orig')
                     self.settings['skip_attr'].append('wake_conv_vel')
 
+
+            if self.settings['save_nonlifting']:
+                self.ClassesToSave += (sharpy.aero.models.nonliftingbodygrid.NonliftingBodyGrid,
+                                       sharpy.utils.datastructures.NonliftingBodyTimeStepInfo,)
+
             if self.settings['save_struct']:
                 self.ClassesToSave += (
                     sharpy.structure.models.beam.Beam,
@@ -222,47 +229,51 @@ def run(self, **kwargs):
 
                     if self.settings['save_struct']:
                         h5utils.add_as_grp(list(),
-                                   hdfile['data']['structure'],
-                                   grpname='timestep_info')
+                                hdfile['data']['structure'],
+                                grpname='timestep_info')
                     if self.settings['save_aero']:
                         h5utils.add_as_grp(list(),
-                                   hdfile['data']['aero'],
-                                   grpname='timestep_info')
+                                hdfile['data']['aero'],
+                                grpname='timestep_info')
+                    if self.settings['save_nonlifting']:
+                        h5utils.add_as_grp(list(),
+                                hdfile['data']['nonlifting_body'],
+                                grpname='timestep_info')
 
                     for it in range(len(self.data.structure.timestep_info)):
                         tstep_p = self.data.structure.timestep_info[it]
                         if tstep_p is not None:
                             self.save_timestep(self.data, self.settings, it, hdfile)
 
-                hdfile.close()
-
-                if self.settings['save_linear_uvlm']:
-                    linhdffile = h5py.File(self.filename.replace('.data.h5', '.uvlmss.h5'), 'a')
-                    h5utils.add_as_grp(self.data.linear.linear_system.uvlm.ss, linhdffile, grpname='ss',
-                                       ClassesToSave=self.ClassesToSave, SkipAttr=self.settings['skip_attr'],
-                                       compress_float=self.settings['compress_float'])
-                    h5utils.add_as_grp(self.data.linear.linear_system.linearisation_vectors, linhdffile,
-                                       grpname='linearisation_vectors',
-                                       ClassesToSave=self.ClassesToSave, SkipAttr=self.settings['skip_attr'],
-                                       compress_float=self.settings['compress_float'])
-                    linhdffile.close()
-
-                if self.settings['save_linear']:
-                    with h5py.File(self.filename_linear, 'a') as linfile:
-                        h5utils.add_as_grp(self.data.linear.linear_system.linearisation_vectors, linfile,
-                                           grpname='linearisation_vectors',
-                                           ClassesToSave=self.ClassesToSave, SkipAttr=self.settings['skip_attr'],
-                                           compress_float=self.settings['compress_float'])
-                        h5utils.add_as_grp(self.data.linear.ss, linfile, grpname='ss',
-                                           ClassesToSave=self.ClassesToSave, SkipAttr=self.settings['skip_attr'],
-                                           compress_float=self.settings['compress_float'])
-
-                if self.settings['save_rom']:
-                    try:
-                        for k, rom in self.data.linear.linear_system.uvlm.rom.items():
-                            rom.save(self.filename.replace('.data.h5', '_{:s}.rom.h5'.format(k.lower())))
-                    except AttributeError:
-                        cout.cout_wrap('Could not locate a reduced order model to save')
+                    hdfile.close()
+
+                    if self.settings['save_linear_uvlm']:
+                        linhdffile = h5py.File(self.filename.replace('.data.h5', '.uvlmss.h5'), 'a')
+                        h5utils.add_as_grp(self.data.linear.linear_system.uvlm.ss, linhdffile, grpname='ss',
+                                        ClassesToSave=self.ClassesToSave, SkipAttr=self.settings['skip_attr'],
+                                        compress_float=self.settings['compress_float'])
+                        h5utils.add_as_grp(self.data.linear.linear_system.linearisation_vectors, linhdffile,
+                                        grpname='linearisation_vectors',
+                                        ClassesToSave=self.ClassesToSave, SkipAttr=self.settings['skip_attr'],
+                                        compress_float=self.settings['compress_float'])
+                        linhdffile.close()
+
+                    if self.settings['save_linear']:
+                        with h5py.File(self.filename_linear, 'a') as linfile:
+                            h5utils.add_as_grp(self.data.linear.linear_system.linearisation_vectors, linfile,
+                                            grpname='linearisation_vectors',
+                                            ClassesToSave=self.ClassesToSave, SkipAttr=self.settings['skip_attr'],
+                                            compress_float=self.settings['compress_float'])
+                            h5utils.add_as_grp(self.data.linear.ss, linfile, grpname='ss',
+                                            ClassesToSave=self.ClassesToSave, SkipAttr=self.settings['skip_attr'],
+                                            compress_float=self.settings['compress_float'])
+
+                    if self.settings['save_rom']:
+                        try:
+                            for k, rom in self.data.linear.linear_system.uvlm.rom.items():
+                                rom.save(self.filename.replace('.data.h5', '_{:s}.rom.h5'.format(k.lower())))
+                        except AttributeError:
+                            cout.cout_wrap('Could not locate a reduced order model to save')
 
             elif self.settings['format'] == 'mat':
                 from scipy.io import savemat
@@ -309,6 +320,13 @@ def save_timestep(data, settings, ts, hdfile):
                                ClassesToSave=(sharpy.utils.datastructures.AeroTimeStepInfo,),
                                SkipAttr=settings['skip_attr'],
                                compress_float=settings['compress_float'])
+        if settings['save_nonlifting']:
+            h5utils.add_as_grp(data.nonlifting_body.timestep_info[ts],
+                               hdfile['data']['nonlifting_body']['timestep_info'],
+                               grpname=("%05d" % ts),
+                               ClassesToSave=(sharpy.utils.datastructures.NonliftingBodyTimeStepInfo,),
+                               SkipAttr=settings['skip_attr'],
+                               compress_float=settings['compress_float'])
         if settings['save_struct']:
             tstep = data.structure.timestep_info[ts]
 
diff --git a/sharpy/postproc/stallcheck.py b/sharpy/postproc/stallcheck.py
index e75fbff29..690cfe43e 100644
--- a/sharpy/postproc/stallcheck.py
+++ b/sharpy/postproc/stallcheck.py
@@ -103,7 +103,7 @@ def check_stall(self):
 
         for i_elem in range(self.data.structure.num_elem):
             for i_local_node in range(self.data.structure.num_node_elem):
-                airfoil_id = self.data.aero.aero_dict['airfoil_distribution'][i_elem, i_local_node]
+                airfoil_id = self.data.aero.data_dict['airfoil_distribution'][i_elem, i_local_node]
                 if self.settings['airfoil_stall_angles']:
                     i_global_node = self.data.structure.connectivities[i_elem, i_local_node]
                     for i_dict in self.data.aero.struct2aero_mapping[i_global_node]:
diff --git a/sharpy/postproc/writevariablestime.py b/sharpy/postproc/writevariablestime.py
index fe532d8df..2a3cd4a8d 100644
--- a/sharpy/postproc/writevariablestime.py
+++ b/sharpy/postproc/writevariablestime.py
@@ -48,6 +48,22 @@ class WriteVariablesTime(BaseSolver):
     settings_default['structure_nodes'] = np.array([-1])
     settings_description['structure_nodes'] = 'Number of the nodes to be writen'
 
+    settings_types['nonlifting_nodes_variables'] = 'list(str)'
+    settings_default['nonlifting_nodes_variables'] = ['']
+    settings_description['nonlifting_nodes_variables'] = 'Variables of :class:`~sharpy.utils.datastructures.NonliftingBodyTimeStepInfo` associated to panels to be writen'
+
+    settings_types['nonlifting_nodes_im'] = 'list(int)'
+    settings_default['nonlifting_nodes_im'] = np.array([0])
+    settings_description['nonlifting_nodes_im'] = 'Chordwise index of the nonlifting panels to be output'
+
+    settings_types['nonlifting_nodes_in'] = 'list(int)'
+    settings_default['nonlifting_nodes_in'] = np.array([0])
+    settings_description['nonlifting_nodes_in'] = 'Spanwise index of the nonlifting panels to be output'
+
+    settings_types['nonlifting_nodes_isurf'] = 'list(int)'
+    settings_default['nonlifting_nodes_isurf'] = np.array([0])
+    settings_description['nonlifting_nodes_isurf'] = "Number of the panels' surface to be output"
+
     settings_types['aero_panels_variables'] = 'list(str)'
     settings_default['aero_panels_variables'] = ['']
     settings_description['aero_panels_variables'] = 'Variables of :class:`~sharpy.utils.datastructures.AeroTimeStepInfo` associated to panels to be writen'
@@ -149,6 +165,17 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
                     if os.path.isfile(filename):
                         os.remove(filename)
 
+        # Nonlifting variables at panels
+        for ivariable in range(len(self.settings['nonlifting_nodes_variables'])):
+            for ipanel in range(len(self.settings['nonlifting_nodes_isurf'])):
+                i_surf = self.settings['nonlifting_nodes_isurf'][ipanel]
+                i_m = self.settings['nonlifting_nodes_im'][ipanel]
+                i_n = self.settings['nonlifting_nodes_in'][ipanel]
+                filename = self.folder + "nonlifting_" + self.settings['nonlifting_nodes_variables'][ivariable] + "_panel" + "_isurf" + str(i_surf) + "_im"+ str(i_m) + "_in"+ str(i_n) + ".dat"
+                if self.settings['cleanup_old_solution']:
+                    if os.path.isfile(filename):
+                        os.remove(filename)
+
         # Aerodynamic variables at panels
         for ivariable in range(len(self.settings['aero_panels_variables'])):
             for ipanel in range(len(self.settings['aero_panels_isurf'])):
@@ -258,6 +285,21 @@ def write(self, it):
                             self.write_nparray_to_file(fid, self.data.ts, var[ielem,inode_in_elem,:], self.settings['delimiter'])
 
 
+        # Aerodynamic variables at nonlifting panels
+        for ivariable in range(len(self.settings['nonlifting_nodes_variables'])):
+            if self.settings['nonlifting_nodes_variables'][ivariable] == '':
+                continue
+            
+            for ipanel in range(len(self.settings['nonlifting_nodes_isurf'])):
+                i_surf = self.settings['nonlifting_nodes_isurf'][ipanel]
+                i_m = self.settings['nonlifting_nodes_im'][ipanel]
+                i_n = self.settings['nonlifting_nodes_in'][ipanel]
+                filename = self.folder + "nonlifting_" + self.settings['nonlifting_nodes_variables'][ivariable] + "_panel" + "_isurf" + str(i_surf) + "_im"+ str(i_m) + "_in"+ str(i_n) + ".dat"
+
+                with open(filename, 'a') as fid:
+                    var = getattr(self.data.nonlifting_body.timestep_info[it], self.settings['nonlifting_nodes_variables'][ivariable])
+                    self.write_value_to_file(fid, self.data.ts, var[i_surf][i_m,i_n], self.settings['delimiter'])
+
         # Aerodynamic variables at panels
         for ivariable in range(len(self.settings['aero_panels_variables'])):
             if self.settings['aero_panels_variables'][ivariable] == '':
@@ -271,7 +313,7 @@ def write(self, it):
 
                 with open(filename, 'a') as fid:
                     var = getattr(self.data.aero.timestep_info[it], self.settings['aero_panels_variables'][ivariable])
-                    self.write_value_to_file(fid, self.data.ts, var.gamma[i_surf][i_m,i_n], self.settings['delimiter'])
+                    self.write_value_to_file(fid, self.data.ts, var[i_surf][i_m,i_n], self.settings['delimiter'])
 
 
         # Aerodynamic variables at nodes
diff --git a/sharpy/solvers/aerogridloader.py b/sharpy/solvers/aerogridloader.py
index bc8f9a8df..af8098676 100644
--- a/sharpy/solvers/aerogridloader.py
+++ b/sharpy/solvers/aerogridloader.py
@@ -6,12 +6,13 @@
 import sharpy.utils.settings as settings_utils
 import sharpy.utils.h5utils as h5utils
 import sharpy.utils.generator_interface as gen_interface
+from sharpy.solvers.gridloader import GridLoader
 
 
 @solver
-class AerogridLoader(BaseSolver):
+class AerogridLoader(GridLoader):
     """
-    ``AerogridLoader`` class, inherited from ``BaseSolver``
+    ``AerogridLoader`` class, inherited from ``GridLoader``
 
     Generates aerodynamic grid based on the input data
 
@@ -36,9 +37,9 @@ class AerogridLoader(BaseSolver):
         settings_types (dict): Acceptable types for the values in ``settings``
         settings_default (dict): Name-value pair of default values for the aerodynamic settings
         data (ProblemData): class structure
-        aero_file_name (str): name of the ``.aero.h5`` HDF5 file
+        file_name (str): name of the ``.aero.h5`` HDF5 file
         aero: empty attribute
-        aero_data_dict (dict): key-value pairs of aerodynamic data
+        data_dict (dict): key-value pairs of aerodynamic data
         wake_shape_generator (class): Wake shape generator
 
     """
@@ -89,28 +90,13 @@ class AerogridLoader(BaseSolver):
                                        settings_options=settings_options)
 
     def __init__(self):
-        self.data = None
-        self.settings = None
-        self.aero_file_name = ''
-        # storage of file contents
-        self.aero_data_dict = dict()
-
-        # aero storage
+        super().__init__
+        self.file_name = '.aero.h5'
         self.aero = None
-
         self.wake_shape_generator = None
 
     def initialise(self, data, restart=False):
-        self.data = data
-        self.settings = data.settings[self.solver_id]
-
-        # init settings
-        settings_utils.to_custom_types(self.settings,
-                           self.settings_types,
-                           self.settings_default, options=self.settings_options)
-
-        # read input file (aero)
-        self.read_files()
+        super().initialise(data)
 
         wake_shape_generator_type = gen_interface.generator_from_string(
             self.settings['wake_shape_generator'])
@@ -119,25 +105,9 @@ def initialise(self, data, restart=False):
                                              self.settings['wake_shape_generator_input'],
                                              restart=restart)
 
-    def read_files(self):
-        # open aero file
-        # first, file names
-        self.aero_file_name = (self.data.case_route +
-                               '/' +
-                               self.data.case_name +
-                               '.aero.h5')
-
-        #  then check that the file exists
-        h5utils.check_file_exists(self.aero_file_name)
-
-        #  read and store the hdf5 file
-        with h5.File(self.aero_file_name, 'r') as aero_file_handle:
-            # store files in dictionary
-            self.aero_data_dict = h5utils.load_h5_in_dict(aero_file_handle)
-
     def run(self, **kwargs):
         self.data.aero = aerogrid.Aerogrid()
-        self.data.aero.generate(self.aero_data_dict,
+        self.data.aero.generate(self.data_dict,
                                 self.data.structure,
                                 self.settings,
                                 self.data.ts)
diff --git a/sharpy/solvers/dynamiccoupled.py b/sharpy/solvers/dynamiccoupled.py
index 72d107467..719a58380 100644
--- a/sharpy/solvers/dynamiccoupled.py
+++ b/sharpy/solvers/dynamiccoupled.py
@@ -178,6 +178,10 @@ class DynamicCoupled(BaseSolver):
                                                  'The dictionary values are dictionaries with the settings ' \
                                                  'needed by each generator.'
 
+    settings_types['nonlifting_body_interactions'] = 'bool'
+    settings_default['nonlifting_body_interactions'] = False
+    settings_description['nonlifting_body_interactions'] = 'Effect of Nonlifting Bodies on Lifting bodies are considered'
+    
     settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
 
@@ -365,17 +369,20 @@ def initialise(self, data, custom_settings=None, restart=False):
 
     def cleanup_timestep_info(self):
         if max(len(self.data.aero.timestep_info), len(self.data.structure.timestep_info)) > 1:
-            # copy last info to first
-            self.data.aero.timestep_info[0] = self.data.aero.timestep_info[-1]
-            self.data.structure.timestep_info[0] = self.data.structure.timestep_info[-1]
-            # delete all the rest
-            while len(self.data.aero.timestep_info) - 1:
-                del self.data.aero.timestep_info[-1]
-            while len(self.data.structure.timestep_info) - 1:
-                del self.data.structure.timestep_info[-1]
+            self.remove_old_timestep_info(self.data.structure.timestep_info)
+            self.remove_old_timestep_info(self.data.aero.timestep_info)            
+            if self.settings['nonlifting_body_interactions']:
+                self.remove_old_timestep_info(self.data.nonlifting_body.timestep_info)
 
         self.data.ts = 0
 
+    def remove_old_timestep_info(self, tstep_info):
+        # copy last info to first
+        tstep_info[0] = tstep_info[-1].copy()
+        # delete all the rest
+        while len(tstep_info) - 1:
+            del tstep_info[-1]
+
     def process_controller_output(self, controlled_state):
         """
         This function modified the solver properties and parameters as
@@ -525,6 +532,10 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None, solvers=No
 
             structural_kstep = self.data.structure.timestep_info[-1].copy()
             aero_kstep = self.data.aero.timestep_info[-1].copy()
+            if self.settings['nonlifting_body_interactions']:
+                nl_body_kstep = self.data.nonlifting_body.timestep_info[-1].copy()
+            else:
+                nl_body_kstep = None
             self.logger.debug('Time step {}'.format(self.data.ts))
 
             # Add the controller here
@@ -565,14 +576,17 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None, solvers=No
                     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)
+                        aero_kstep,
+                        nl_body_kstep)
                     break
 
                 # generate new grid (already rotated)
                 aero_kstep = controlled_aero_kstep.copy()
+                
                 self.aero_solver.update_custom_grid(
-                    structural_kstep,
-                    aero_kstep)
+                        structural_kstep,
+                        aero_kstep,
+                        nl_body_kstep)
 
                 # compute unsteady contribution
                 force_coeff = 0.0
@@ -604,7 +618,8 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None, solvers=No
                 self.data = self.aero_solver.run(aero_step=aero_kstep,
                                                  structural_step=structural_kstep,
                                                  convect_wake=True,
-                                                 unsteady_contribution=unsteady_contribution)
+                                                 unsteady_contribution=unsteady_contribution,
+                                                 nl_body_tstep = nl_body_kstep)
                 self.time_aero += time.perf_counter() - ini_time_aero
 
                 previous_kstep = structural_kstep.copy()
@@ -615,12 +630,15 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None, solvers=No
                 previous_kstep.runtime_unsteady_forces = previous_runtime_unsteady_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.aero_solver.update_custom_grid(
+                        structural_kstep,
+                        aero_kstep,
+                        nl_body_kstep)
+                
                 self.map_forces(aero_kstep,
-                                structural_kstep,
-                                force_coeff)
+                            structural_kstep,
+                            nl_body_kstep = nl_body_kstep,
+                            unsteady_forces_coeff = force_coeff)
 
                 # relaxation
                 relax_factor = self.relaxation_factor(k)
@@ -664,16 +682,20 @@ def time_loop(self, in_queue=None, out_queue=None, finish_event=None, solvers=No
                                     self.aero_solver,
                                     self.with_runtime_generators):
                     # move the aerodynamic surface according to the structural one
-                    self.aero_solver.update_custom_grid(
-                        structural_kstep,
-                        aero_kstep)
+                    self.aero_solver.update_custom_grid(structural_kstep,
+                                                        aero_kstep,
+                                                        nl_body_tstep = nl_body_kstep)
                     break
 
             # move the aerodynamic surface according the the structural one
-            self.aero_solver.update_custom_grid(structural_kstep, aero_kstep)
+            self.aero_solver.update_custom_grid(structural_kstep,
+                                                aero_kstep,
+                                                nl_body_tstep = nl_body_kstep)
 
             self.aero_solver.add_step()
             self.data.aero.timestep_info[-1] = aero_kstep.copy()
+            if self.settings['nonlifting_body_interactions']:
+                self.data.nonlifting_body.timestep_info[-1] = nl_body_kstep.copy()
             self.structural_solver.add_step()
             self.data.structure.timestep_info[-1] = structural_kstep.copy()
 
@@ -748,7 +770,9 @@ def convergence(self, k, tstep, previous_tstep,
             return False
 
         # Check the special case of no aero and no runtime generators
-        if aero_solver.solver_id.lower() == "noaero" and not with_runtime_generators:
+        if (aero_solver.solver_id.lower() == "noaero"\
+             or struct_solver.solver_id.lower()  == "nostructural")\
+            and not with_runtime_generators:
             return True
 
         # relative residuals
@@ -787,7 +811,7 @@ def convergence(self, k, tstep, previous_tstep,
                         return True
 
 
-    def map_forces(self, aero_kstep, structural_kstep, unsteady_forces_coeff=1.0):
+    def map_forces(self, aero_kstep, structural_kstep, nl_body_kstep = None, unsteady_forces_coeff=1.0):
         # set all forces to 0
         structural_kstep.steady_applied_forces.fill(0.0)
         structural_kstep.unsteady_applied_forces.fill(0.0)
@@ -802,8 +826,8 @@ def map_forces(self, aero_kstep, structural_kstep, unsteady_forces_coeff=1.0):
             self.data.structure.node_master_elem,
             self.data.structure.connectivities,
             structural_kstep.cag(),
-            self.data.aero.aero_dict)
-        dynamic_struct_forces = mapping.aero2struct_force_mapping(
+            self.data.aero.data_dict)
+        dynamic_struct_forces = unsteady_forces_coeff*mapping.aero2struct_force_mapping(
             aero_kstep.dynamic_forces,
             self.data.aero.struct2aero_mapping,
             aero_kstep.zeta,
@@ -812,7 +836,7 @@ def map_forces(self, aero_kstep, structural_kstep, unsteady_forces_coeff=1.0):
             self.data.structure.node_master_elem,
             self.data.structure.connectivities,
             structural_kstep.cag(),
-            self.data.aero.aero_dict)
+            self.data.aero.data_dict)
 
         if self.correct_forces:
             struct_forces = \
@@ -825,6 +849,18 @@ def map_forces(self, aero_kstep, structural_kstep, unsteady_forces_coeff=1.0):
         structural_kstep.postproc_node['aero_steady_forces'] = struct_forces
         structural_kstep.postproc_node['aero_unsteady_forces'] = dynamic_struct_forces
 
+        # if self.settings['nonlifting_body_interactions']:
+        #     struct_forces +=  mapping.aero2struct_force_mapping(
+        #         nl_body_kstep.forces,
+        #         self.data.nonlifting_body.struct2aero_mapping,
+        #         nl_body_kstep.zeta,
+        #         structural_kstep.pos,
+        #         structural_kstep.psi,
+        #         self.data.structure.node_master_elem,
+        #         self.data.structure.connectivities,
+        #         structural_kstep.cag(),
+        #         self.data.nonlifting_body.data_dict)
+        # prescribed forces + aero forces
         # prescribed forces + aero forces + runtime generated
         structural_kstep.steady_applied_forces += struct_forces
         structural_kstep.steady_applied_forces += self.data.structure.ini_info.steady_applied_forces
diff --git a/sharpy/solvers/gridloader.py b/sharpy/solvers/gridloader.py
new file mode 100644
index 000000000..0caad5dc4
--- /dev/null
+++ b/sharpy/solvers/gridloader.py
@@ -0,0 +1,64 @@
+import h5py as h5
+import numpy as np
+
+from sharpy.utils.solver_interface import solver, BaseSolver
+import sharpy.utils.settings as settings_utils
+import sharpy.utils.h5utils as h5utils
+
+
+@solver
+class GridLoader(BaseSolver):
+    """
+    ``GridLoader`` class, inherited from ``BaseSolver``
+
+    Parent class for Aerogridloader and Nonliftingbodygridloader. Both classes
+    generate aerodynamic grids based on the input data
+
+    Args:
+        data (PreSharpy): ``ProblemData`` class structure
+
+    Attributes:
+        settings (dict): Name-value pair of the settings employed by the aerodynamic solver
+        settings_types (dict): Acceptable types for the values in ``settings``
+        settings_default (dict): Name-value pair of default values for the aerodynamic settings
+        data (ProblemData): class structure
+        afile_name (str): name of the HDF5 file, e.g. ``.aero.h5``
+        aero: empty attribute
+        data_dict (dict): key-value pairs of aerodynamic data
+
+    """
+    solver_id = 'GridLoader'
+    solver_classification = 'other'
+
+    settings_types = dict()
+    settings_default = dict()
+    settings_description = dict()
+    settings_options = dict()
+
+    def __init__(self):
+        self.data = None
+        self.settings = None
+        self.file_name = ''
+        self.data_dict = dict()
+
+    def initialise(self, data, restart=False):
+        self.data = data
+        self.read_input_files()
+        
+        self.settings = data.settings[self.solver_id]
+        settings_utils.to_custom_types(self.settings,
+                                       self.settings_types,
+                                       self.settings_default, options=self.settings_options)
+
+
+    def read_input_files(self):
+        self.file_name = (self.data.case_route +
+                          '/' +
+                          self.data.case_name +
+                          self.file_name)
+                                                    
+        h5utils.check_file_exists(self.file_name)
+
+        #  read and store the hdf5 file in dictionary
+        with h5.File(self.file_name, 'r') as file_handle:
+            self.data_dict = h5utils.load_h5_in_dict(file_handle)
\ No newline at end of file
diff --git a/sharpy/solvers/nonliftingbodygridloader.py b/sharpy/solvers/nonliftingbodygridloader.py
new file mode 100644
index 000000000..85f80bac2
--- /dev/null
+++ b/sharpy/solvers/nonliftingbodygridloader.py
@@ -0,0 +1,44 @@
+from sharpy.utils.solver_interface import solver
+import sharpy.aero.models.nonliftingbodygrid as nonliftingbodygrid
+import sharpy.utils.settings as settings_utils
+from sharpy.solvers.gridloader import GridLoader
+
+
+@solver
+class NonliftingbodygridLoader(GridLoader):
+    """
+    ``NonliftingbodygridLoader`` class, inherited from ``GridLoader``
+
+    Generates aerodynamic grid for nonlifting bodies based on the input data
+
+    Args:
+        data (PreSharpy): ``ProblemData`` class structure
+
+    Attributes:
+        settings (dict): Name-value pair of the settings employed by the aerodynamic solver
+        settings_types (dict): Acceptable types for the values in ``settings``
+        settings_default (dict): Name-value pair of default values for the aerodynamic settings
+        data (ProblemData): class structure
+        file_name (str): name of the ``.nonlifting_body.h5`` HDF5 file
+        aero: empty attribute
+        aero_data_dict (dict): key-value pairs of aerodynamic data
+
+
+    """
+    solver_id = 'NonliftingbodygridLoader'
+    solver_classification = 'loader'
+
+    def __init__(self):
+        super().__init__
+        self.file_name = '.nonlifting_body.h5'
+
+        # nonlifting_body storage
+        self.nonlifting_body = None
+
+    def run(self, **kwargs):
+        self.data.nonlifting_body = nonliftingbodygrid.NonliftingBodyGrid()
+        self.data.nonlifting_body.generate(self.data_dict,
+                                           self.data.structure,
+                                           self.settings,
+                                           self.data.ts)
+        return self.data
diff --git a/sharpy/solvers/nonlineardynamiccoupledstep.py b/sharpy/solvers/nonlineardynamiccoupledstep.py
index 621686ba5..c6a665b83 100644
--- a/sharpy/solvers/nonlineardynamiccoupledstep.py
+++ b/sharpy/solvers/nonlineardynamiccoupledstep.py
@@ -77,8 +77,7 @@ def initialise(self, data, custom_settings=None, restart=False):
     def run(self, **kwargs):
 
         structural_step = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
-        # TODO: previous_structural_step never used
-        previous_structural_step = settings_utils.set_value_or_default(kwargs, 'previous_structural_step', self.data.structure.timestep_info[-1])
+        
         dt= settings_utils.set_value_or_default(kwargs, 'dt', self.settings['dt'])
 
         xbeamlib.xbeam_step_couplednlndyn(self.data.structure,
diff --git a/sharpy/solvers/staticcoupled.py b/sharpy/solvers/staticcoupled.py
index 8fcb8a343..e3fde1bae 100644
--- a/sharpy/solvers/staticcoupled.py
+++ b/sharpy/solvers/staticcoupled.py
@@ -74,6 +74,10 @@ class StaticCoupled(BaseSolver):
     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_types['nonlifting_body_interactions'] = 'bool'
+    settings_default['nonlifting_body_interactions'] = False
+    settings_description['nonlifting_body_interactions'] = 'Consider forces induced by nonlifting bodies'
 
     settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
@@ -149,17 +153,20 @@ def increase_ts(self):
 
     def cleanup_timestep_info(self):
         if max(len(self.data.aero.timestep_info), len(self.data.structure.timestep_info)) > 1:
-            # copy last info to first
-            self.data.aero.timestep_info[0] = self.data.aero.timestep_info[-1].copy()
-            self.data.structure.timestep_info[0] = self.data.structure.timestep_info[-1].copy()
-            # delete all the rest
-            while len(self.data.aero.timestep_info) - 1:
-                del self.data.aero.timestep_info[-1]
-            while len(self.data.structure.timestep_info) - 1:
-                del self.data.structure.timestep_info[-1]
+            self.remove_old_timestep_info(self.data.structure.timestep_info)
+            self.remove_old_timestep_info(self.data.aero.timestep_info)            
+            if self.settings['nonlifting_body_interactions']:
+                self.remove_old_timestep_info(self.data.nonlifting_body.timestep_info)
 
         self.data.ts = 0
 
+    def remove_old_timestep_info(self, tstep_info):
+        # copy last info to first
+        tstep_info[0] = tstep_info[-1].copy()
+        # delete all the rest
+        while len(tstep_info) - 1:
+            del tstep_info[-1]
+
     def run(self, **kwargs):
         for i_step in range(self.settings['n_load_steps'] + 1):
             if (i_step == self.settings['n_load_steps'] and
@@ -189,14 +196,29 @@ def run(self, **kwargs):
                     self.data.structure.node_master_elem,
                     self.data.structure.connectivities,
                     self.data.structure.timestep_info[self.data.ts].cag(),
-                    self.data.aero.aero_dict)
-
+                    self.data.aero.data_dict)
+                        
                 if self.correct_forces:
                     struct_forces = \
                         self.correct_forces_generator.generate(aero_kstep=self.data.aero.timestep_info[self.data.ts],
                                                                structural_kstep=self.data.structure.timestep_info[self.data.ts],
                                                                struct_forces=struct_forces,
                                                                ts=0)
+
+                # map nonlifting forces to structural nodes
+                if self.settings['nonlifting_body_interactions']:
+                    struct_forces += mapping.aero2struct_force_mapping(
+                        self.data.nonlifting_body.timestep_info[self.data.ts].forces,
+                        self.data.nonlifting_body.struct2aero_mapping,
+                        self.data.nonlifting_body.timestep_info[self.data.ts].zeta,
+                        self.data.structure.timestep_info[self.data.ts].pos,
+                        self.data.structure.timestep_info[self.data.ts].psi,
+                        self.data.structure.node_master_elem,
+                        self.data.structure.connectivities,
+                        self.data.structure.timestep_info[self.data.ts].cag(),
+                        self.data.nonlifting_body.data_dict,
+                        skip_moments_generated_by_forces = True)
+
                 self.data.aero.timestep_info[self.data.ts].aero_steady_forces_beam_dof = struct_forces
                 self.data.structure.timestep_info[self.data.ts].postproc_node['aero_steady_forces'] = struct_forces  # B
                 
@@ -239,6 +261,7 @@ def run(self, **kwargs):
                 # update grid
                 self.aero_solver.update_step()
 
+                self.structural_solver.update(self.data.structure.timestep_info[self.data.ts])
                 # convergence
                 if self.convergence(i_iter, i_step):
                     # create q and dqdt vectors
@@ -342,7 +365,7 @@ def change_trim(self, alpha, thrust, thrust_nodes, tail_deflection, tail_cs_inde
 
         # tail deflection
         try:
-            self.data.aero.aero_dict['control_surface_deflection'][tail_cs_index] = tail_deflection
+            self.data.aero.data_dict['control_surface_deflection'][tail_cs_index] = tail_deflection
         except KeyError:
             raise Exception('This model has no control surfaces')
         except IndexError:
diff --git a/sharpy/solvers/staticuvlm.py b/sharpy/solvers/staticuvlm.py
index 66fcd7d35..e9fa9b859 100644
--- a/sharpy/solvers/staticuvlm.py
+++ b/sharpy/solvers/staticuvlm.py
@@ -43,6 +43,18 @@ class StaticUvlm(BaseSolver):
     settings_default['horseshoe'] = False
     settings_description['horseshoe'] = 'Horseshoe wake modelling for steady simulations.'
 
+    settings_types['nonlifting_body_interactions'] = 'bool'
+    settings_default['nonlifting_body_interactions'] = False
+    settings_description['nonlifting_body_interactions'] = 'Consider nonlifting body interactions'
+
+    settings_types['only_nonlifting'] = 'bool'
+    settings_default['only_nonlifting'] = False
+    settings_description['only_nonlifting'] = 'Consider only nonlifting bodies'
+
+    settings_types['phantom_wing_test'] = 'bool'
+    settings_default['phantom_wing_test'] = False
+    settings_description['phantom_wing_test'] = 'Debug option'
+
     settings_types['num_cores'] = 'int'
     settings_default['num_cores'] = 0
     settings_description['num_cores'] = 'Number of cores to use in the VLM lib'
@@ -111,6 +123,10 @@ class StaticUvlm(BaseSolver):
     settings_default['map_forces_on_struct'] = False
     settings_description['map_forces_on_struct'] = 'Maps the forces on the structure at the end of the timestep. Only usefull if the solver is used outside StaticCoupled'
 
+    settings_types['ignore_first_x_nodes_in_force_calculation'] = 'int'
+    settings_default['ignore_first_x_nodes_in_force_calculation'] = 0
+    settings_description['ignore_first_x_nodes_in_force_calculation'] = 'Ignores the forces on the first user-specified number of nodes of all surfaces.'
+
     settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description)
 
@@ -138,60 +154,84 @@ def initialise(self, data, custom_settings=None, restart=False):
 
     def add_step(self):
         self.data.aero.add_timestep()
+        if self.settings['nonlifting_body_interactions']:
+            self.data.nonlifting_body.add_timestep()
+            
 
     def update_grid(self, beam):
-        self.data.aero.generate_zeta(beam,
-                                     self.data.aero.aero_settings,
-                                     -1,
-                                     beam_ts=-1)
 
-    def update_custom_grid(self, structure_tstep, aero_tstep):
+        if not self.settings['only_nonlifting']:
+            self.data.aero.generate_zeta(beam,
+                                        self.data.aero.aero_settings,
+                                        -1,
+                                        beam_ts=-1)
+        if self.settings['nonlifting_body_interactions'] or self.settings['only_nonlifting']:
+            self.data.nonlifting_body.generate_zeta(beam,
+                                                    self.data.nonlifting_body.aero_settings,
+                                                    -1,
+                                                    beam_ts=-1)
+
+    def update_custom_grid(self, structure_tstep, aero_tstep, nonlifting_tstep=None):
         self.data.aero.generate_zeta_timestep_info(structure_tstep,
                                                    aero_tstep,
                                                    self.data.structure,
                                                    self.data.aero.aero_settings,
                                                    dt=self.settings['rollup_dt'])
+        if self.settings['nonlifting_body_interactions']:
+            self.data.nonlifting_body.generate_zeta_timestep_info(structure_tstep,
+                                                    nonlifting_tstep,
+                                                    self.data.structure,
+                                                    self.data.nonlifting_body.aero_settings)
 
     def run(self, **kwargs):
 
-        aero_tstep = settings_utils.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[-1])
-        structure_tstep = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[-1])
-        dt = settings_utils.set_value_or_default(kwargs, 'dt', self.settings['rollup_dt'])
-        t = settings_utils.set_value_or_default(kwargs, 't', self.data.ts*dt)
-
-        unsteady_contribution = False
-        convect_wake = False
-
-        if not aero_tstep.zeta:
-            return self.data
-
-        # generate the wake because the solid shape might change
-        self.data.aero.wake_shape_generator.generate({'zeta': aero_tstep.zeta,
-                                            'zeta_star': aero_tstep.zeta_star,
-                                            'gamma': aero_tstep.gamma,
-                                            'gamma_star': aero_tstep.gamma_star,
-                                            'dist_to_orig': aero_tstep.dist_to_orig})
-
-        # generate uext
-        self.velocity_generator.generate({'zeta': aero_tstep.zeta,
-                                          'override': True,
-                                          'for_pos': structure_tstep.for_pos[0:3]},
-                                          aero_tstep.u_ext)
-        # grid orientation
-        uvlmlib.vlm_solver(aero_tstep,
-                           self.settings)
-
-        if self.settings['map_forces_on_struct']:
-            structure_tstep.steady_applied_forces[:] = mapping.aero2struct_force_mapping(
-                    aero_tstep.forces,
-                    self.data.aero.struct2aero_mapping,
-                    self.data.aero.timestep_info[self.data.ts].zeta,
-                    structure_tstep.pos,
-                    structure_tstep.psi,
-                    self.data.structure.node_master_elem,
-                    self.data.structure.connectivities,
-                    structure_tstep.cag(),
-                    self.data.aero.aero_dict)
+        structure_tstep = settings_utils.set_value_or_default(kwargs, 'structural_step', self.data.structure.timestep_info[self.data.ts])
+        
+        if not self.settings['only_nonlifting']:
+            aero_tstep = settings_utils.set_value_or_default(kwargs, 'aero_step', self.data.aero.timestep_info[self.data.ts])
+            if not self.data.aero.timestep_info[self.data.ts].zeta:
+                return self.data
+
+            # generate the wake because the solid shape might change
+            self.data.aero.wake_shape_generator.generate({'zeta': aero_tstep.zeta,
+                                                'zeta_star': aero_tstep.zeta_star,
+                                                'gamma': aero_tstep.gamma,
+                                                'gamma_star': aero_tstep.gamma_star,
+                                                'dist_to_orig': aero_tstep.dist_to_orig})
+        
+            if self.settings['nonlifting_body_interactions']:
+                # generate uext
+                self.velocity_generator.generate({'zeta': self.data.nonlifting_body.timestep_info[self.data.ts].zeta,
+                                                'override': True,
+                                                'for_pos': structure_tstep.for_pos[0:3]},
+                                                self.data.nonlifting_body.timestep_info[self.data.ts].u_ext)
+                # generate uext
+                self.velocity_generator.generate({'zeta': self.data.aero.timestep_info[self.data.ts].zeta,
+                                                'override': True,
+                                                'for_pos': self.data.structure.timestep_info[self.data.ts].for_pos[0:3]},  
+                                                self.data.aero.timestep_info[self.data.ts].u_ext)
+                # grid orientation
+                uvlmlib.vlm_solver_lifting_and_nonlifting_bodies(self.data.aero.timestep_info[self.data.ts],
+                                                            self.data.nonlifting_body.timestep_info[self.data.ts],
+                                                            self.settings)
+            else:        
+                # generate uext
+                self.velocity_generator.generate({'zeta': self.data.aero.timestep_info[self.data.ts].zeta,
+                                                'override': True,
+                                                        'for_pos': self.data.structure.timestep_info[self.data.ts].for_pos[0:3]},
+                                                        self.data.aero.timestep_info[self.data.ts].u_ext)
+
+
+                # grid orientation
+                uvlmlib.vlm_solver(self.data.aero.timestep_info[self.data.ts],
+                                    self.settings)
+        else:
+            self.velocity_generator.generate({'zeta': self.data.nonlifting_body.timestep_info[self.data.ts].zeta,
+                                            'override': True,
+                                            'for_pos': self.data.structure.timestep_info[self.data.ts].for_pos[0:3]},
+                                            self.data.nonlifting_body.timestep_info[self.data.ts].u_ext)
+            uvlmlib.vlm_solver_nonlifting_body(self.data.nonlifting_body.timestep_info[self.data.ts],
+                                            self.settings)
 
         return self.data
 
@@ -199,12 +239,17 @@ def next_step(self):
         """ Updates de aerogrid based on the info of the step, and increases
         the self.ts counter """
         self.data.aero.add_timestep()
+        if self.settings['nonlifting_body_interactions']:
+            self.data.nonlifting_body.add_timestep()
         self.update_step()
 
     def update_step(self):
-        self.data.aero.generate_zeta(self.data.structure,
-                                     self.data.aero.aero_settings,
-                                     self.data.ts)
-        # for i_surf in range(self.data.aero.timestep_info[self.data.ts].n_surf):
-        #     self.data.aero.timestep_info[self.data.ts].forces[i_surf].fill(0.0)
-        #     self.data.aero.timestep_info[self.data.ts].dynamic_forces[i_surf].fill(0.0)
+        if not self.settings['only_nonlifting']:
+            self.data.aero.generate_zeta(self.data.structure,
+                                        self.data.aero.aero_settings,
+                                        self.data.ts)
+        if self.settings['nonlifting_body_interactions'] or self.settings['only_nonlifting']:
+            self.data.nonlifting_body.generate_zeta(self.data.structure,
+                                                    self.data.nonlifting_body.aero_settings,
+                                                    self.data.ts)
+
diff --git a/sharpy/solvers/stepuvlm.py b/sharpy/solvers/stepuvlm.py
index d4cd9ad40..ba4b2a686 100644
--- a/sharpy/solvers/stepuvlm.py
+++ b/sharpy/solvers/stepuvlm.py
@@ -130,6 +130,26 @@ class StepUvlm(BaseSolver):
     settings_types['quasi_steady'] = 'bool'
     settings_default['quasi_steady'] = False
     settings_description['quasi_steady'] = 'Use quasi-steady approximation in UVLM'
+    
+    settings_types['only_nonlifting'] = 'bool'
+    settings_default['only_nonlifting'] = False
+    settings_description['only_nonlifting'] = 'Consider nonlifting body interactions'
+
+    settings_types['nonlifting_body_interactions'] = 'bool'
+    settings_default['nonlifting_body_interactions'] = False
+    settings_description['nonlifting_body_interactions'] = 'Consider nonlifting body interactions'
+
+    settings_types['phantom_wing_test'] = 'bool'
+    settings_default['phantom_wing_test'] = False
+    settings_description['phantom_wing_test'] = 'Debug option'
+
+    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_types['ignore_first_x_nodes_in_force_calculation'] = 'int'
+    settings_default['ignore_first_x_nodes_in_force_calculation'] = 0
+    settings_description['ignore_first_x_nodes_in_force_calculation'] = 'Ignores the forces on the first user-specified number of nodes of all surfaces.'
 
     settings_table = settings_utils.SettingsTable()
     __doc__ += settings_table.generate(settings_types, settings_default, settings_description, settings_options)
@@ -218,13 +238,32 @@ def run(self, **kwargs):
                                               'for_pos': structure_tstep.for_pos,
                                               'is_wake': True},
                                              aero_tstep.u_ext_star)
+        if self.settings['nonlifting_body_interactions']:
 
-        uvlmlib.uvlm_solver(self.data.ts,
-                            aero_tstep,
-                            structure_tstep,
-                            self.settings,
-                            convect_wake=convect_wake,
-                            dt=dt)
+            nl_body_tstep = settings_utils.set_value_or_default(kwargs, 'nl_body_tstep', self.data.nonlifting_body.timestep_info[-1])
+            self.velocity_generator.generate({'zeta': nl_body_tstep.zeta,
+                                              'override': True,
+                                              'ts': self.data.ts,
+                                              'dt': dt,
+                                              't': t,
+                                              'for_pos': structure_tstep.for_pos,
+                                              'is_wake': False},
+                                             nl_body_tstep.u_ext)
+            
+            uvlmlib.uvlm_solver_lifting_and_nonlifting(self.data.ts,
+                                                       aero_tstep, 
+                                                       nl_body_tstep,
+                                                       structure_tstep, 
+                                                       self.settings,
+                                                       convect_wake=convect_wake, 
+                                                       dt=dt)
+        else:            
+            uvlmlib.uvlm_solver(self.data.ts,
+                                aero_tstep,
+                                structure_tstep,
+                                self.settings,
+                                convect_wake=convect_wake,
+                                dt=dt)
 
         if unsteady_contribution and not self.settings['quasi_steady']:
             # calculate unsteady (added mass) forces:
@@ -248,24 +287,38 @@ def run(self, **kwargs):
         else:
             for i_surf in range(len(aero_tstep.gamma)):
                 aero_tstep.gamma_dot[i_surf][:] = 0.0
-
         return self.data
 
     def add_step(self):
         self.data.aero.add_timestep()
+        if self.settings['nonlifting_body_interactions']:              
+            self.data.nonlifting_body.add_timestep()
 
     def update_grid(self, beam):
         self.data.aero.generate_zeta(beam,
                                      self.data.aero.aero_settings,
                                      -1,
                                      beam_ts=-1)
+        if self.settings['nonlifting_body_interactions']:                            
+            self.data.nonlifting_body.generate_zeta(beam,
+                                                    self.data.aero.aero_settings,
+                                                    -1,
+                                                    beam_ts=-1)
 
-    def update_custom_grid(self, structure_tstep, aero_tstep):
+    def update_custom_grid(self, structure_tstep, aero_tstep, nl_body_tstep = None):
         self.data.aero.generate_zeta_timestep_info(structure_tstep,
                                                    aero_tstep,
                                                    self.data.structure,
                                                    self.data.aero.aero_settings,
                                                    dt=self.settings['dt'])
+        if self.settings['nonlifting_body_interactions']:
+            if nl_body_tstep is None:
+                nl_body_tstep  =  self.data.nonlifting_body.timestep_info[-1] 
+            self.data.nonlifting_body.generate_zeta_timestep_info(structure_tstep,
+                                                                  nl_body_tstep,
+                                                                  self.data.structure,
+                                                                  self.data.nonlifting_body.aero_settings,
+                                                                  dt = self.settings['dt'])
 
     @staticmethod
     def filter_gamma_dot(tstep, history, filter_param):
diff --git a/sharpy/solvers/trim.py b/sharpy/solvers/trim.py
index b8543f21e..4cf0a5a05 100644
--- a/sharpy/solvers/trim.py
+++ b/sharpy/solvers/trim.py
@@ -291,7 +291,7 @@ def solver_wrapper(x, x_info, solver_data, i_dim=-1):
     tstep.quat[:] = orientation_quat
     # control surface deflection
     for i_cs in range(len(x_info['i_control_surfaces'])):
-        solver_data.data.aero.aero_dict['control_surface_deflection'][x_info['control_surfaces_id'][i_cs]] = x[x_info['i_control_surfaces'][i_cs]]
+        solver_data.data.aero.data_dict['control_surface_deflection'][x_info['control_surfaces_id'][i_cs]] = x[x_info['i_control_surfaces'][i_cs]]
     # thrust input
     tstep.steady_applied_forces[:] = 0.0
     try:
diff --git a/sharpy/structure/utils/modalutils.py b/sharpy/structure/utils/modalutils.py
index bdcbd7781..f16353c26 100644
--- a/sharpy/structure/utils/modalutils.py
+++ b/sharpy/structure/utils/modalutils.py
@@ -166,8 +166,8 @@ def get_mode_zeta(data, eigvect):
             # print('%.2d,%.2d'%(nn,ss))
 
             # surface panelling
-            M = aero.aero_dimensions[ss][0]
-            N = aero.aero_dimensions[ss][1]
+            M = aero.dimensions[ss][0]
+            N = aero.dimensions[ss][1]
 
             for mm in range(M + 1):
                 # get position of vertex in B FoR
diff --git a/sharpy/utils/datastructures.py b/sharpy/utils/datastructures.py
index 57b0156fe..29222ad3d 100644
--- a/sharpy/utils/datastructures.py
+++ b/sharpy/utils/datastructures.py
@@ -11,21 +11,19 @@
 import sharpy.utils.multibody as mb
 
 
-class AeroTimeStepInfo(object):
+class TimeStepInfo(object):
     """
-    Aerodynamic Time step class.
+    Time step class.
 
-    Contains the relevant aerodynamic attributes for a single time step. All variables should be expressed in ``G``
-    FoR unless otherwise stated.
+    It is the parent class of the AeroTimeStepInfo and NonliftingBodyTimeStepInfo, which contain the relevant
+    aerodynamic attributes for a single time step. All variables should be expressed in ``G`` FoR unless
+    otherwise stated.
 
     Attributes:
         ct_dimensions: Pointer to ``dimensions`` to interface the C++ library `uvlmlib``
-        ct_dimensions_star: Pointer to ``dimensions_star`` to interface the C++ library `uvlmlib``
 
         dimensions (np.ndarray): Matrix defining the dimensions of the vortex grid on solid surfaces
           ``[num_surf x chordwise panels x spanwise panels]``
-        dimensions_star (np.ndarray): Matrix defining the dimensions of the vortex grid on wakes
-          ``[num_surf x streamwise panels x spanwise panels]``
 
         n_surf (int): Number of aerodynamic surfaces on solid bodies. Each aerodynamic surface on solid bodies will
           have an associted wake.
@@ -39,17 +37,8 @@ class AeroTimeStepInfo(object):
           ``[n_surf][3 x (chordwise nodes + 1) x (spanwise nodes + 1)]``
         dynamic_forces (list(np.ndarray)): Forces associated to time derivatives on grid vertices
           ``[n_surf][3 x (chordwise nodes + 1) x (spanwise nodes + 1)]``
-        zeta_star (list(np.ndarray): Location of wake grid vertices
-          ``[n_surf][3 x (streamwise nodes + 1) x (spanwise nodes + 1)]``
         u_ext (list(np.ndarray)): Background flow velocity on solid grid nodes
           ``[n_surf][3 x (chordwise nodes + 1) x (spanwise nodes + 1)]``
-        u_ext_star (list(np.ndarray)): Background flow velocity on wake grid nodes
-          ``[n_surf][3 x (streamwise nodes + 1) x (spanwise nodes + 1)]``
-        gamma (list(np.ndarray)): Circulation associated to solid panels
-          ``[n_surf][3 x chordwise nodes x spanwise nodes]``
-        gamma_star (list(np.ndarray)): Circulation associated to wake panels
-          ``[n_surf][3 x streamwise nodes x spanwise nodes]``
-        gamma_dot (list(np.ndarray)): Time derivative of ``gamma``
 
         inertial_steady_forces (list(np.ndarray)): Total aerodynamic steady forces in ``G`` FoR ``[n_surf x 6]``
         body_steady_forces (list(np.ndarray)): Total aerodynamic steady forces in ``A`` FoR ``[n_surf x 6]``
@@ -62,20 +51,15 @@ class AeroTimeStepInfo(object):
         in_global_AFoR (bool): ``True`` if the variables are stored in the global A FoR. ``False`` if they are stored
           in the local A FoR of each body. Always ``True`` for single-body simulations. Currently not used.
 
-        control_surface_deflection (np.ndarray): Deflection of the control surfaces, in `rad` and if fitted.
 
     Args:
         dimensions (np.ndarray): Matrix defining the dimensions of the vortex grid on solid surfaces
           ``[num_surf x chordwise panels x spanwise panels]``
-        dimensions_star (np.ndarray): Matrix defining the dimensions of the vortex grid on wakes
-          ``[num_surf x streamwise panels x spanwise panels]``
     """
-    def __init__(self, dimensions, dimensions_star):
+    def __init__(self, dimensions):
         self.ct_dimensions = None
-        self.ct_dimensions_star = None
 
         self.dimensions = dimensions.copy()
-        self.dimensions_star = dimensions_star.copy()
         self.n_surf = self.dimensions.shape[0]
         # generate placeholder for aero grid zeta coordinates
         self.zeta = []
@@ -95,10 +79,9 @@ def __init__(self, dimensions, dimensions_star):
         self.normals = []
         for i_surf in range(self.n_surf):
             self.normals.append(np.zeros((3,
-                                          dimensions[i_surf, 0],
-                                          dimensions[i_surf, 1]),
-                                         dtype=ct.c_double))
-
+                                         dimensions[i_surf, 0],
+                                         dimensions[i_surf, 1]),
+                                        dtype=ct.c_double))
         # panel forces
         self.forces = []
         for i_surf in range(self.n_surf):
@@ -114,14 +97,6 @@ def __init__(self, dimensions, dimensions_star):
                                                  dimensions[i_surf, 1] + 1),
                                                 dtype=ct.c_double))
 
-        # generate placeholder for aero grid zeta_star coordinates
-        self.zeta_star = []
-        for i_surf in range(self.n_surf):
-            self.zeta_star.append(np.zeros((3,
-                                            dimensions_star[i_surf, 0] + 1,
-                                            dimensions_star[i_surf, 1] + 1),
-                                           dtype=ct.c_double))
-
         # placeholder for external velocity
         self.u_ext = []
         for i_surf in range(self.n_surf):
@@ -130,39 +105,6 @@ def __init__(self, dimensions, dimensions_star):
                                         dimensions[i_surf, 1] + 1),
                                        dtype=ct.c_double))
 
-        self.u_ext_star = []
-        for i_surf in range(self.n_surf):
-            self.u_ext_star.append(np.zeros((3,
-                                             dimensions_star[i_surf, 0] + 1,
-                                             dimensions_star[i_surf, 1] + 1),
-                                            dtype=ct.c_double))
-
-        # allocate gamma and gamma star matrices
-        self.gamma = []
-        for i_surf in range(self.n_surf):
-            self.gamma.append(np.zeros((dimensions[i_surf, 0],
-                                        dimensions[i_surf, 1]),
-                                       dtype=ct.c_double))
-
-        self.gamma_star = []
-        for i_surf in range(self.n_surf):
-            self.gamma_star.append(np.zeros((dimensions_star[i_surf, 0],
-                                             dimensions_star[i_surf, 1]),
-                                            dtype=ct.c_double))
-
-        self.gamma_dot = []
-        for i_surf in range(self.n_surf):
-            self.gamma_dot.append(np.zeros((dimensions[i_surf, 0],
-                                            dimensions[i_surf, 1]),
-                                           dtype=ct.c_double))
-
-        # Distance from the trailing edge of the wake vertices
-        self.dist_to_orig = []
-        for i_surf in range(self.n_surf):
-            self.dist_to_orig.append(np.zeros((dimensions_star[i_surf, 0] + 1,
-                                               dimensions_star[i_surf, 1] + 1),
-                                               dtype=ct.c_double))
-
         # total forces - written by AeroForcesCalculator
         self.inertial_steady_forces = np.zeros((self.n_surf, 6))
         self.body_steady_forces = np.zeros((self.n_surf, 6))
@@ -175,13 +117,14 @@ def __init__(self, dimensions, dimensions_star):
         # Multibody variables
         self.in_global_AFoR = True
 
-        self.control_surface_deflection = np.array([])
 
     def copy(self):
         """
-        Returns a copy of a deepcopy of a :class:`~sharpy.utils.datastructures.AeroTimeStepInfo`
+        Returns a copy of a deepcopy of a :class:`~sharpy.utils.datastructures.TimeStepInfo`
         """
-        copied = AeroTimeStepInfo(self.dimensions, self.dimensions_star)
+        return self.create_placeholder(TimeStepInfo(self.dimensions))
+
+    def create_placeholder(self, copied):
         # generate placeholder for aero grid zeta coordinates
         for i_surf in range(copied.n_surf):
             copied.zeta[i_surf] = self.zeta[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
@@ -201,30 +144,10 @@ def copy(self):
         for i_surf in range(copied.n_surf):
             copied.dynamic_forces[i_surf] = self.dynamic_forces[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
 
-        # generate placeholder for aero grid zeta_star coordinates
-        for i_surf in range(copied.n_surf):
-            copied.zeta_star[i_surf] = self.zeta_star[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
-
         # placeholder for external velocity
         for i_surf in range(copied.n_surf):
             copied.u_ext[i_surf] = self.u_ext[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
 
-        for i_surf in range(copied.n_surf):
-            copied.u_ext_star[i_surf] = self.u_ext_star[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
-
-        # allocate gamma and gamma star matrices
-        for i_surf in range(copied.n_surf):
-            copied.gamma[i_surf] = self.gamma[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
-
-        for i_surf in range(copied.n_surf):
-            copied.gamma_dot[i_surf] = self.gamma_dot[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
-
-        for i_surf in range(copied.n_surf):
-            copied.gamma_star[i_surf] = self.gamma_star[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
-
-        for i_surf in range(copied.n_surf):
-            copied.dist_to_orig[i_surf] = self.dist_to_orig[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
-
         # total forces
         copied.inertial_steady_forces = self.inertial_steady_forces.astype(dtype=ct.c_double, copy=True, order='C')
         copied.body_steady_forces = self.body_steady_forces.astype(dtype=ct.c_double, copy=True, order='C')
@@ -234,8 +157,6 @@ def copy(self):
         copied.postproc_cell = copy.deepcopy(self.postproc_cell)
         copied.postproc_node = copy.deepcopy(self.postproc_node)
 
-        copied.control_surface_deflection = self.control_surface_deflection.astype(dtype=ct.c_double, copy=True)
-
         return copied
 
     def generate_ctypes_pointers(self):
@@ -243,7 +164,6 @@ def generate_ctypes_pointers(self):
         Generates the pointers to aerodynamic variables used to interface the C++ library ``uvlmlib``
         """
         self.ct_dimensions = self.dimensions.astype(dtype=ct.c_uint, copy=True)
-        self.ct_dimensions_star = self.dimensions_star.astype(dtype=ct.c_uint, copy=True)
 
         n_surf = len(self.dimensions)
 
@@ -259,33 +179,11 @@ def generate_ctypes_pointers(self):
             for i_dim in range(NDIM):
                 self.ct_zeta_dot_list.append(self.zeta_dot[i_surf][i_dim, :, :].reshape(-1))
 
-        self.ct_zeta_star_list = []
-        for i_surf in range(self.n_surf):
-            for i_dim in range(NDIM):
-                self.ct_zeta_star_list.append(self.zeta_star[i_surf][i_dim, :, :].reshape(-1))
-
         self.ct_u_ext_list = []
         for i_surf in range(self.n_surf):
             for i_dim in range(NDIM):
                 self.ct_u_ext_list.append(self.u_ext[i_surf][i_dim, :, :].reshape(-1))
 
-        self.ct_u_ext_star_list = []
-        for i_surf in range(self.n_surf):
-            for i_dim in range(NDIM):
-                self.ct_u_ext_star_list.append(self.u_ext_star[i_surf][i_dim, :, :].reshape(-1))
-
-        self.ct_gamma_list = []
-        for i_surf in range(self.n_surf):
-            self.ct_gamma_list.append(self.gamma[i_surf][:, :].reshape(-1))
-
-        self.ct_gamma_dot_list = []
-        for i_surf in range(self.n_surf):
-            self.ct_gamma_dot_list.append(self.gamma_dot[i_surf][:, :].reshape(-1))
-
-        self.ct_gamma_star_list = []
-        for i_surf in range(self.n_surf):
-            self.ct_gamma_star_list.append(self.gamma_star[i_surf][:, :].reshape(-1))
-
         self.ct_normals_list = []
         for i_surf in range(self.n_surf):
             for i_dim in range(NDIM):
@@ -301,10 +199,6 @@ def generate_ctypes_pointers(self):
             for i_dim in range(NDIM*2):
                 self.ct_dynamic_forces_list.append(self.dynamic_forces[i_surf][i_dim, :, :].reshape(-1))
 
-        self.ct_dist_to_orig_list = []
-        for i_surf in range(self.n_surf):
-            self.ct_dist_to_orig_list.append(self.dist_to_orig[i_surf][:, :].reshape(-1))
-
         try:
             self.postproc_cell['incidence_angle']
         except KeyError:
@@ -319,32 +213,18 @@ def generate_ctypes_pointers(self):
 
         self.ct_p_dimensions = ((ct.POINTER(ct.c_uint)*n_surf)
                                 (* np.ctypeslib.as_ctypes(self.ct_dimensions)))
-        self.ct_p_dimensions_star = ((ct.POINTER(ct.c_uint)*n_surf)
-                                     (* np.ctypeslib.as_ctypes(self.ct_dimensions_star)))
         self.ct_p_zeta = ((ct.POINTER(ct.c_double)*len(self.ct_zeta_list))
                           (* [np.ctypeslib.as_ctypes(array) for array in self.ct_zeta_list]))
         self.ct_p_zeta_dot = ((ct.POINTER(ct.c_double)*len(self.ct_zeta_dot_list))
                           (* [np.ctypeslib.as_ctypes(array) for array in self.ct_zeta_dot_list]))
-        self.ct_p_zeta_star = ((ct.POINTER(ct.c_double)*len(self.ct_zeta_star_list))
-                               (* [np.ctypeslib.as_ctypes(array) for array in self.ct_zeta_star_list]))
         self.ct_p_u_ext = ((ct.POINTER(ct.c_double)*len(self.ct_u_ext_list))
                            (* [np.ctypeslib.as_ctypes(array) for array in self.ct_u_ext_list]))
-        self.ct_p_u_ext_star = ((ct.POINTER(ct.c_double)*len(self.ct_u_ext_star_list))
-                           (* [np.ctypeslib.as_ctypes(array) for array in self.ct_u_ext_star_list]))
-        self.ct_p_gamma = ((ct.POINTER(ct.c_double)*len(self.ct_gamma_list))
-                           (* [np.ctypeslib.as_ctypes(array) for array in self.ct_gamma_list]))
-        self.ct_p_gamma_dot = ((ct.POINTER(ct.c_double)*len(self.ct_gamma_dot_list))
-                               (* [np.ctypeslib.as_ctypes(array) for array in self.ct_gamma_dot_list]))
-        self.ct_p_gamma_star = ((ct.POINTER(ct.c_double)*len(self.ct_gamma_star_list))
-                                (* [np.ctypeslib.as_ctypes(array) for array in self.ct_gamma_star_list]))
         self.ct_p_normals = ((ct.POINTER(ct.c_double)*len(self.ct_normals_list))
                              (* [np.ctypeslib.as_ctypes(array) for array in self.ct_normals_list]))
         self.ct_p_forces = ((ct.POINTER(ct.c_double)*len(self.ct_forces_list))
                             (* [np.ctypeslib.as_ctypes(array) for array in self.ct_forces_list]))
         self.ct_p_dynamic_forces = ((ct.POINTER(ct.c_double)*len(self.ct_dynamic_forces_list))
                             (* [np.ctypeslib.as_ctypes(array) for array in self.ct_dynamic_forces_list]))
-        self.ct_p_dist_to_orig = ((ct.POINTER(ct.c_double)*len(self.ct_dist_to_orig_list))
-                           (* [np.ctypeslib.as_ctypes(array) for array in self.ct_dist_to_orig_list]))
 
         if with_incidence_angle:
             self.postproc_cell['incidence_angle_ct_pointer'] = ((ct.POINTER(ct.c_double)*len(self.ct_incidence_list))
@@ -352,83 +232,357 @@ def generate_ctypes_pointers(self):
 
     def remove_ctypes_pointers(self):
         """
-        Removes the pointers to aerodynamic variables used to interface the C++ library ``uvlmlib``
+            Removes the pointers to aerodynamic variables used to interface the C++ library ``uvlmlib``
         """
-        try:
-            del self.ct_p_dimensions
-        except AttributeError:
-            pass
 
-        try:
-            del self.ct_p_dimensions_star
-        except AttributeError:
-            pass
+        list_class_attributes = list(self.__dict__.keys()).copy()
+        for name_attribute in list_class_attributes:
+            if "ct_p_" in name_attribute:
+                self.__delattr__(name_attribute)
 
-        try:
-            del self.ct_p_zeta
-        except AttributeError:
-            pass
+        for k in list(self.postproc_cell.keys()):
+            if 'ct_list' in k:
+                del self.postproc_cell[k]
+            elif 'ct_pointer' in k:
+                del self.postproc_cell[k]
 
-        try:
-            del self.ct_p_zeta_star
-        except AttributeError:
-            pass
 
-        try:
-            del self.ct_p_zeta_dot
-        except AttributeError:
-            pass
+class NonliftingBodyTimeStepInfo(TimeStepInfo):
+    """
+    Nonlifting Body Time step class.
 
-        try:
-            del self.ct_p_u_ext
-        except AttributeError:
-            pass
+    It is the inheritance from ``TimeStepInfo`` and contains the relevant aerodynamic attributes for
+    a single time step of a nonlifting body. All variables should be expressed in ``G``
+    FoR unless otherwise stated.
 
-        try:
-            del self.ct_p_u_ext_star
-        except AttributeError:
-            pass
+    Attributes:
+        ct_dimensions: Pointer to ``dimensions`` to interface the C++ library `uvlmlib``
 
-        try:
-            del self.ct_p_gamma
-        except AttributeError:
-            pass
+        dimensions (np.ndarray): Matrix defining the dimensions of the vortex grid on solid surfaces
+        ``[num_surf x radial panels x spanwise panels]``
 
-        try:
-            del self.ct_p_gamma_dot
-        except AttributeError:
-            pass
+        n_surf (int): Number of aerodynamic surfaces on nonlifting bodies.
 
-        try:
-            del self.ct_p_gamma_star
-        except AttributeError:
-            pass
+        zeta (list(np.ndarray): Location of solid grid vertices
+          ``[n_surf][3 x (radial panel) x (spanwise panel)]``
+        zeta_dot (list(np.ndarray)): Time derivative of ``zeta``
+        normals (list(np.ndarray)): Normal direction to panels at the panel center
+          ``[n_surf][3 x radial panels x spanwise panels]``
+        forces (list(np.ndarray)): Forces not associated to time derivatives on grid vertices
+          ``[n_surf][3 x (radial panels) x (spanwise panels)]``
+        dynamic_forces (list(np.ndarray)): Forces associated to time derivatives on grid vertices
+          ``[n_surf][3 x (radial panels) x (spanwise panels)]``
+
+        u_ext (list(np.ndarray)): Background flow velocity on solid grid panel
+          ``[n_surf][3 x (radial panels) x (spanwise panel + 1)]``
+        sigma (list(np.ndarray)): Source strength associated to solid panels
+          ``[n_surf][3 x radial panel x spanwise panel]``
+        sigma_dot (list(np.ndarray)): Time derivative of ``sigma``
+        pressure_coefficients (list(np.ndarray)): Pressure coefficient associated to solid panels
+          ``[n_surf][radial panel x spanwise panel]``
+
+        inertial_total_forces (list(np.ndarray)): Total aerodynamic forces in ``G`` FoR ``[n_surf x 6]``
+        body_total_forces (list(np.ndarray)): Total aerodynamic forces in ``A`` FoR ``[n_surf x 6]``
+        inertial_steady_forces (list(np.ndarray)): Total aerodynamic steady forces in ``G`` FoR ``[n_surf x 6]``
+        body_steady_forces (list(np.ndarray)): Total aerodynamic steady forces in ``A`` FoR ``[n_surf x 6]``
+        inertial_unsteady_forces (list(np.ndarray)): Total aerodynamic unsteady forces in ``G`` FoR ``[n_surf x 6]``
+        body_unsteady_forces (list(np.ndarray)): Total aerodynamic unsteady forces in ``A`` FoR ``[n_surf x 6]``
 
-        try:
-            del self.ct_p_normals
-        except AttributeError:
-            pass
+        postproc_cell (dict): Variables associated to cells to be postprocessed
+        postproc_node (dict): Variables associated to panel to be postprocessed
 
-        try:
-            del self.ct_p_forces
-        except AttributeError:
-            pass
+        in_global_AFoR (bool): ``True`` if the variables are stored in the global A FoR. ``False`` if they are stored
+          in the local A FoR of each body. Always ``True`` for single-body simulations. Currently not used.
 
-        try:
-            del self.ct_p_dynamic_forces
-        except AttributeError:
-            pass
+    Args:
+        dimensions (np.ndarray): Matrix defining the dimensions of the vortex grid on solid surfaces
+          ``[num_surf x radial panels x spanwise panels]``
+    """
+    def __init__(self, dimensions): #remove dimensions_star as input
+        super().__init__(dimensions)
 
-        try:
-            del self.ct_p_dist_to_orig
-        except AttributeError:
-            pass
+        # allocate sigma matrices
+        self.sigma = []
+        for i_surf in range(self.n_surf):
+            self.sigma.append(np.zeros((dimensions[i_surf, 0],
+                                        dimensions[i_surf, 1]),
+                                       dtype=ct.c_double))
 
-        for k in list(self.postproc_cell.keys()):
-            if 'ct_list' in k:
-                del self.postproc_cell[k]
-            elif 'ct_pointer' in k:
-                del self.postproc_cell[k]
+        self.sigma_dot = []
+        for i_surf in range(self.n_surf):
+            self.sigma_dot.append(np.zeros((dimensions[i_surf, 0],
+                                            dimensions[i_surf, 1]),
+                                           dtype=ct.c_double))
+
+        self.pressure_coefficients = []
+        for i_surf in range(self.n_surf):
+            self.pressure_coefficients.append(np.zeros((dimensions[i_surf, 0],
+                                        dimensions[i_surf, 1]),
+                                        dtype=ct.c_double))
+    def copy(self):
+        """
+        Returns a copy of a deepcopy of a :class:`~sharpy.utils.datastructures.AeroTimeStepInfo`
+        """
+        return self.create_placeholder(NonliftingBodyTimeStepInfo(self.dimensions))
+
+    def create_placeholder(self, copied):
+        super().create_placeholder(copied)
+
+        # allocate sigma matrices
+        for i_surf in range(copied.n_surf):
+            copied.sigma[i_surf] = self.sigma[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+        for i_surf in range(copied.n_surf):
+            copied.sigma_dot[i_surf] = self.sigma_dot[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+        for i_surf in range(copied.n_surf):
+            copied.pressure_coefficients[i_surf] = self.pressure_coefficients[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+
+        return copied
+
+
+    def generate_ctypes_pointers(self):
+        """
+        Generates the pointers to aerodynamic variables used to interface the C++ library ``uvlmlib``
+        """
+        super().generate_ctypes_pointers()
+
+        from sharpy.utils.constants import NDIM
+
+        self.ct_sigma_list = []
+        for i_surf in range(self.n_surf):
+            self.ct_sigma_list.append(self.sigma[i_surf][:, :].reshape(-1))
+
+        self.ct_sigma_dot_list = []
+        for i_surf in range(self.n_surf):
+            self.ct_sigma_dot_list.append(self.sigma_dot[i_surf][:, :].reshape(-1))
+
+        self.ct_pressure_coefficients_list = []
+        for i_surf in range(self.n_surf):
+            self.ct_pressure_coefficients_list.append(self.pressure_coefficients[i_surf][:, :].reshape(-1))
+
+        self.ct_p_sigma = ((ct.POINTER(ct.c_double)*len(self.ct_sigma_list))
+                            (* [np.ctypeslib.as_ctypes(array) for array in self.ct_sigma_list]))
+        self.ct_p_sigma_dot = ((ct.POINTER(ct.c_double)*len(self.ct_sigma_dot_list))
+                            (* [np.ctypeslib.as_ctypes(array) for array in self.ct_sigma_list]))
+        self.ct_p_pressure_coefficients = ((ct.POINTER(ct.c_double)*len(self.ct_pressure_coefficients_list))
+                                (* [np.ctypeslib.as_ctypes(array) for array in self.ct_pressure_coefficients_list]))
+
+
+
+class AeroTimeStepInfo(TimeStepInfo):
+    """
+    Aerodynamic Time step class.
+
+    Contains the relevant aerodynamic attributes for a single time step. All variables should be expressed in ``G``
+    FoR unless otherwise stated.
+
+    Attributes:
+        ct_dimensions: Pointer to ``dimensions`` to interface the C++ library `uvlmlib``
+        ct_dimensions_star: Pointer to ``dimensions_star`` to interface the C++ library `uvlmlib``
+
+        dimensions (np.ndarray): Matrix defining the dimensions of the vortex grid on solid surfaces
+          ``[num_surf x chordwise panels x spanwise panels]``
+        dimensions_star (np.ndarray): Matrix defining the dimensions of the vortex grid on wakes
+          ``[num_surf x streamwise panels x spanwise panels]``
+
+        n_surf (int): Number of aerodynamic surfaces on solid bodies. Each aerodynamic surface on solid bodies will
+          have an associted wake.
+
+        zeta (list(np.ndarray): Location of solid grid vertices
+          ``[n_surf][3 x (chordwise nodes + 1) x (spanwise nodes + 1)]``
+        zeta_dot (list(np.ndarray)): Time derivative of ``zeta``
+        normals (list(np.ndarray)): Normal direction to panels at the panel center
+          ``[n_surf][3 x chordwise nodes x spanwise nodes]``
+        forces (list(np.ndarray)): Forces not associated to time derivatives on grid vertices
+          ``[n_surf][3 x (chordwise nodes + 1) x (spanwise nodes + 1)]``
+        dynamic_forces (list(np.ndarray)): Forces associated to time derivatives on grid vertices
+          ``[n_surf][3 x (chordwise nodes + 1) x (spanwise nodes + 1)]``
+        zeta_star (list(np.ndarray): Location of wake grid vertices
+          ``[n_surf][3 x (streamwise nodes + 1) x (spanwise nodes + 1)]``
+        u_ext (list(np.ndarray)): Background flow velocity on solid grid nodes
+          ``[n_surf][3 x (chordwise nodes + 1) x (spanwise nodes + 1)]``
+        u_ext_star (list(np.ndarray)): Background flow velocity on wake grid nodes
+          ``[n_surf][3 x (streamwise nodes + 1) x (spanwise nodes + 1)]``
+        gamma (list(np.ndarray)): Circulation associated to solid panels
+          ``[n_surf][3 x chordwise nodes x spanwise nodes]``
+        gamma_star (list(np.ndarray)): Circulation associated to wake panels
+          ``[n_surf][3 x streamwise nodes x spanwise nodes]``
+        gamma_dot (list(np.ndarray)): Time derivative of ``gamma``
+
+        inertial_total_forces (list(np.ndarray)): Total aerodynamic forces in ``G`` FoR ``[n_surf x 6]``
+        body_total_forces (list(np.ndarray)): Total aerodynamic forces in ``A`` FoR ``[n_surf x 6]``
+        inertial_steady_forces (list(np.ndarray)): Total aerodynamic steady forces in ``G`` FoR ``[n_surf x 6]``
+        body_steady_forces (list(np.ndarray)): Total aerodynamic steady forces in ``A`` FoR ``[n_surf x 6]``
+        inertial_unsteady_forces (list(np.ndarray)): Total aerodynamic unsteady forces in ``G`` FoR ``[n_surf x 6]``
+        body_unsteady_forces (list(np.ndarray)): Total aerodynamic unsteady forces in ``A`` FoR ``[n_surf x 6]``
+
+        postproc_cell (dict): Variables associated to cells to be postprocessed
+        postproc_node (dict): Variables associated to nodes to be postprocessed
+
+        in_global_AFoR (bool): ``True`` if the variables are stored in the global A FoR. ``False`` if they are stored
+          in the local A FoR of each body. Always ``True`` for single-body simulations. Currently not used.
+
+        control_surface_deflection (np.ndarray): Deflection of the control surfaces, in `rad` and if fitted.
+
+    Args:
+        dimensions (np.ndarray): Matrix defining the dimensions of the vortex grid on solid surfaces
+          ``[num_surf x chordwise panels x spanwise panels]``
+        dimensions_star (np.ndarray): Matrix defining the dimensions of the vortex grid on wakes
+          ``[num_surf x streamwise panels x spanwise panels]``
+    """
+    def __init__(self, dimensions, dimensions_star):
+        super().__init__(dimensions)
+        self.ct_dimensions_star = None
+
+        self.dimensions_star = dimensions_star.copy()
+
+        # generate placeholder for aero grid zeta_star coordinates
+        self.zeta_star = []
+        for i_surf in range(self.n_surf):
+            self.zeta_star.append(np.zeros((3,
+                                            dimensions_star[i_surf, 0] + 1,
+                                            dimensions_star[i_surf, 1] + 1),
+                                           dtype=ct.c_double))
+
+        self.u_ext_star = []
+        for i_surf in range(self.n_surf):
+            self.u_ext_star.append(np.zeros((3,
+                                             dimensions_star[i_surf, 0] + 1,
+                                             dimensions_star[i_surf, 1] + 1),
+                                            dtype=ct.c_double))
+
+        # allocate gamma and gamma star matrices
+        self.gamma = []
+        for i_surf in range(self.n_surf):
+            self.gamma.append(np.zeros((dimensions[i_surf, 0],
+                                        dimensions[i_surf, 1]),
+                                       dtype=ct.c_double))
+
+        self.gamma_star = []
+        for i_surf in range(self.n_surf):
+            self.gamma_star.append(np.zeros((dimensions_star[i_surf, 0],
+                                             dimensions_star[i_surf, 1]),
+                                            dtype=ct.c_double))
+
+        self.gamma_dot = []
+        for i_surf in range(self.n_surf):
+            self.gamma_dot.append(np.zeros((dimensions[i_surf, 0],
+                                            dimensions[i_surf, 1]),
+                                           dtype=ct.c_double))
+
+        # Distance from the trailing edge of the wake vertices
+        self.dist_to_orig = []
+        for i_surf in range(self.n_surf):
+            self.dist_to_orig.append(np.zeros((dimensions_star[i_surf, 0] + 1,
+                                               dimensions_star[i_surf, 1] + 1),
+                                               dtype=ct.c_double))
+
+        self.wake_conv_vel = []
+        for i_surf in range(self.n_surf):
+            self.wake_conv_vel.append(np.zeros((dimensions_star[i_surf, 0],
+                                            dimensions_star[i_surf, 1]),
+                                           dtype=ct.c_double))
+
+        # Junction handling
+        self.flag_zeta_phantom = np.zeros((1, self.n_surf),
+                                            dtype=ct.c_int)
+        self.control_surface_deflection = np.array([])
+
+    def copy(self):
+        return self.create_placeholder(AeroTimeStepInfo(self.dimensions, self.dimensions_star))
+
+    def create_placeholder(self, copied):
+        super().create_placeholder(copied)
+        # generate placeholder for aero grid zeta_star coordinates
+        for i_surf in range(copied.n_surf):
+            copied.zeta_star[i_surf] = self.zeta_star[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+        # placeholder for external velocity
+        for i_surf in range(copied.n_surf):
+            copied.u_ext_star[i_surf] = self.u_ext_star[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+        # allocate gamma and gamma star matrices
+        for i_surf in range(copied.n_surf):
+            copied.gamma[i_surf] = self.gamma[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+        for i_surf in range(copied.n_surf):
+            copied.gamma_dot[i_surf] = self.gamma_dot[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+        for i_surf in range(copied.n_surf):
+            copied.gamma_star[i_surf] = self.gamma_star[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+        for i_surf in range(copied.n_surf):
+            copied.dist_to_orig[i_surf] = self.dist_to_orig[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+        for i_surf in range(copied.n_surf):
+            copied.wake_conv_vel[i_surf] = self.wake_conv_vel[i_surf].astype(dtype=ct.c_double, copy=True, order='C')
+
+        copied.control_surface_deflection = self.control_surface_deflection.astype(dtype=ct.c_double, copy=True)
+
+        # phantom panel flags
+        copied.flag_zeta_phantom = self.flag_zeta_phantom.astype(dtype=ct.c_int, copy=True, order='C')
+        
+        return copied
+
+    def generate_ctypes_pointers(self):
+        from sharpy.utils.constants import NDIM
+        n_surf = len(self.dimensions)
+        super().generate_ctypes_pointers()
+        self.ct_dimensions_star = self.dimensions_star.astype(dtype=ct.c_uint, copy=True)
+
+        self.ct_zeta_star_list = []
+        for i_surf in range(self.n_surf):
+            for i_dim in range(NDIM):
+                self.ct_zeta_star_list.append(self.zeta_star[i_surf][i_dim, :, :].reshape(-1))
+
+        self.ct_u_ext_star_list = []
+        for i_surf in range(self.n_surf):
+            for i_dim in range(NDIM):
+                self.ct_u_ext_star_list.append(self.u_ext_star[i_surf][i_dim, :, :].reshape(-1))
+
+        self.ct_gamma_list = []
+        for i_surf in range(self.n_surf):
+            self.ct_gamma_list.append(self.gamma[i_surf][:, :].reshape(-1))
+
+        self.ct_gamma_dot_list = []
+        for i_surf in range(self.n_surf):
+            self.ct_gamma_dot_list.append(self.gamma_dot[i_surf][:, :].reshape(-1))
+
+        self.ct_gamma_star_list = []
+        for i_surf in range(self.n_surf):
+            self.ct_gamma_star_list.append(self.gamma_star[i_surf][:, :].reshape(-1))
+
+        self.ct_dist_to_orig_list = []
+        for i_surf in range(self.n_surf):
+            self.ct_dist_to_orig_list.append(self.dist_to_orig[i_surf][:, :].reshape(-1))
+
+        self.ct_wake_conv_vel_list = []
+        for i_surf in range(self.n_surf):
+            self.ct_wake_conv_vel_list.append(self.wake_conv_vel[i_surf][:, :].reshape(-1))
+        self.ct_flag_zeta_phantom_list = self.flag_zeta_phantom[:].reshape(-1)
+ 
+
+
+        self.ct_p_dimensions_star = ((ct.POINTER(ct.c_uint)*n_surf)
+                                     (* np.ctypeslib.as_ctypes(self.ct_dimensions_star)))
+        self.ct_p_zeta_star = ((ct.POINTER(ct.c_double)*len(self.ct_zeta_star_list))
+                               (* [np.ctypeslib.as_ctypes(array) for array in self.ct_zeta_star_list]))
+        self.ct_p_u_ext_star = ((ct.POINTER(ct.c_double)*len(self.ct_u_ext_star_list))
+                           (* [np.ctypeslib.as_ctypes(array) for array in self.ct_u_ext_star_list]))
+        self.ct_p_gamma = ((ct.POINTER(ct.c_double)*len(self.ct_gamma_list))
+                           (* [np.ctypeslib.as_ctypes(array) for array in self.ct_gamma_list]))
+        self.ct_p_gamma_dot = ((ct.POINTER(ct.c_double)*len(self.ct_gamma_dot_list))
+                               (* [np.ctypeslib.as_ctypes(array) for array in self.ct_gamma_dot_list]))
+        self.ct_p_gamma_star = ((ct.POINTER(ct.c_double)*len(self.ct_gamma_star_list))
+                                (* [np.ctypeslib.as_ctypes(array) for array in self.ct_gamma_star_list]))
+        self.ct_p_dist_to_orig = ((ct.POINTER(ct.c_double)*len(self.ct_dist_to_orig_list))
+                           (* [np.ctypeslib.as_ctypes(array) for array in self.ct_dist_to_orig_list]))
+        self.ct_p_wake_conv_vel = ((ct.POINTER(ct.c_double)*len(self.ct_wake_conv_vel_list))
+                           (* [np.ctypeslib.as_ctypes(array) for array in self.ct_wake_conv_vel_list]))
+        self.ct_p_flag_zeta_phantom = np.ctypeslib.as_ctypes(self.ct_flag_zeta_phantom_list)
+ 
 
 
 def init_matrix_structure(dimensions, with_dim_dimension, added_size=0):
diff --git a/sharpy/utils/generate_cases.py b/sharpy/utils/generate_cases.py
index e2b489fe2..bbcd295f3 100644
--- a/sharpy/utils/generate_cases.py
+++ b/sharpy/utils/generate_cases.py
@@ -1768,7 +1768,9 @@ def set_default_values(self):
             else:
                 solver_name = solver
             self.solvers[solver_name] = deepcopy(aux_solvers[solver])
-            
+            if solver in ['GridLoader', 'NonliftingbodygridLoader']:
+                # Skip this solver as no default values for GridLoader exist. 
+                continue
 
     def check(self):
 
diff --git a/sharpy/utils/generator_interface.py b/sharpy/utils/generator_interface.py
index ff3b827f9..70fd3a28d 100644
--- a/sharpy/utils/generator_interface.py
+++ b/sharpy/utils/generator_interface.py
@@ -28,10 +28,7 @@ def print_available_generators():
 
 
 class BaseGenerator(metaclass=ABCMeta):
-    
-    def teardown(self):
-        pass
-
+    pass
 
 def generator_from_string(string):
     return dict_of_generators[string]
diff --git a/sharpy/utils/h5utils.py b/sharpy/utils/h5utils.py
index 927af014c..5f3549f83 100644
--- a/sharpy/utils/h5utils.py
+++ b/sharpy/utils/h5utils.py
@@ -69,7 +69,7 @@ def check_fem_dict(fem_dict):
         print(' PASSED')
 
 
-def check_aero_dict(aero_dict):
+def check_data_dict(data_dict):
     pass
 
 
diff --git a/sharpy/utils/solver_interface.py b/sharpy/utils/solver_interface.py
index 6703c258d..85dae79eb 100644
--- a/sharpy/utils/solver_interface.py
+++ b/sharpy/utils/solver_interface.py
@@ -103,8 +103,11 @@ def dictionary_of_solvers(print_info=True):
     import sharpy.postproc
     dictionary = dict()
     for solver in dict_of_solvers:
-        init_solver = initialise_solver(solver, print_info)
-        dictionary[solver] = init_solver.settings_default
+        if solver not in ['GridLoader', 'NonliftingBodyGridLoader']:
+            init_solver = initialise_solver(solver, print_info)
+            dictionary[solver] = init_solver.settings_default
+        else:
+            dictionary[solver] = {}
 
     return dictionary
 
diff --git a/tests/coupled/dynamic/hale/generate_hale.py b/tests/coupled/dynamic/hale/generate_hale.py
index b7399f273..2d68873bb 100644
--- a/tests/coupled/dynamic/hale/generate_hale.py
+++ b/tests/coupled/dynamic/hale/generate_hale.py
@@ -10,18 +10,13 @@
 # EXECUTION
 flow = ['BeamLoader',
         'AerogridLoader',
-        'StaticTrim',
+        'StaticCoupled',
         'DynamicCoupled',
         'BeamLoads'
         ]
 
 # if free_flight is False, the motion of the centre of the wing is prescribed.
 free_flight = True
-if not free_flight:
-    case_name += '_prescribed'
-    amplitude = 0 * np.pi / 180
-    period = 3
-    case_name += '_amp_' + str(amplitude).replace('.', '') + '_period_' + str(period)
 
 # FLIGHT CONDITIONS
 # the simulation is set such that the aircraft flies at a u_inf velocity while
@@ -44,7 +39,7 @@
 # gust settings
 gust_intensity = 0.20
 gust_length = 1 * u_inf
-gust_offset = 0.2 * u_inf
+gust_offset = 0.0 * u_inf
 
 # numerics
 n_step = 5
@@ -101,7 +96,7 @@
 # chordiwse panels
 m = 4
 # spanwise elements
-n_elem_multiplier = 2
+n_elem_multiplier = 1
 n_elem_main = int(4 * n_elem_multiplier)
 n_elem_tail = int(2 * n_elem_multiplier)
 n_elem_fin = int(2 * n_elem_multiplier)
@@ -112,7 +107,7 @@
 physical_time = 30
 tstep_factor = 1.
 dt = 1.0 / m / u_inf * tstep_factor
-n_tstep = 20
+n_tstep = 5
 
 # END OF INPUT-----------------------------------------------------------------
 
@@ -669,15 +664,16 @@ def generate_solver_file():
                                   'structural_solver_settings': settings[solver],
                                   'aero_solver': 'StepUvlm',
                                   'aero_solver_settings': settings['StepUvlm'],
-                                  'fsi_substeps': 200,
-                                  'fsi_tolerance': fsi_tolerance,
-                                  'relaxation_factor': relaxation_factor,
+                                  'fsi_substeps': 3,
+                                  'fsi_tolerance': 1e-3,
+                                  'relaxation_factor': 0,
                                   'minimum_steps': 1,
                                   'relaxation_steps': 150,
                                   'final_relaxation_factor': 0.5,
                                   'n_time_steps': n_tstep,
                                   'dt': dt,
-                                  'include_unsteady_force_contribution': 'on',}
+                                  'include_unsteady_force_contribution': 'on',
+                                  }
     
     settings['BeamLoader'] = {'unsteady': 'on',
                               'orientation': algebra.euler2quat(np.array([roll,
diff --git a/tests/coupled/dynamic/test_dynamic.py b/tests/coupled/dynamic/test_dynamic.py
index 99cbca807..635832c6f 100644
--- a/tests/coupled/dynamic/test_dynamic.py
+++ b/tests/coupled/dynamic/test_dynamic.py
@@ -1,8 +1,6 @@
 import numpy as np
-import importlib
 import unittest
 import os
-import sharpy.utils.cout_utils as cout
 
 
 class TestCoupledDynamic(unittest.TestCase):
@@ -12,33 +10,31 @@ class TestCoupledDynamic(unittest.TestCase):
     - Gust response of the hale aircraft
     """
 
-    @classmethod
-    def setUpClass(cls):
-        # run all the cases generators
-        case = 'hale'
-        mod = importlib.import_module('tests.coupled.dynamic.' + case + '.generate_' + case)
-        pass
-
+    route_file_dir = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))
+    
     def test_hale_dynamic(self):
         """
         Case and results from:
         tests/coupled/dynamic/hale 
-        reference results produced with SHARPy version 1.3
+        reference results produced with SHARPy version 2.0
         :return:
         """
         import sharpy.sharpy_main
-
+        try:
+            import hale.generate_hale
+        except:
+            import tests.coupled.dynamic.hale.generate_hale
+        
         case_name = 'hale'
-        route_test_dir = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))
-        cases_folder = os.path.join(route_test_dir, case_name)
+        cases_folder = os.path.join(self.route_file_dir, case_name)
         output_folder = cases_folder + '/output/'
 
         sharpy.sharpy_main.main(['', cases_folder + '/hale.sharpy'])
-        n_tstep = 20
-
+        n_tstep = 5
         # compare results with reference values 
-        ref_Fz = 50.4986064826483
-        ref_My = -1833.91402522644
+        ref_Fz = -531.023900359779
+        ref_My = -1530.0477841197576
+
         file = os.path.join(output_folder, case_name, 'beam/beam_loads_%i.csv' % (n_tstep))
         beam_loads_ts = np.loadtxt(file, delimiter=',')
         np.testing.assert_almost_equal(float(beam_loads_ts[0, 6]), ref_Fz,
@@ -51,18 +47,21 @@ def test_hale_dynamic(self):
                                        verbose=True)
 
     @classmethod
-    def tearDownClass(cls):
-
+    def tearDown(self):
+        """
+            Removes all created files within this test.
+        """
         import shutil
-        list_cases = ['hale']
-        list_file_extensions = ['.fem.h5', '.aero.h5', '.sharpy']
-        list_folders = ['output', '__pycache__']
-        for case in list_cases:
-            file_path = os.path.join(os.path.abspath(os.path.dirname(os.path.realpath(__file__))),
-                                     case)
-            for folder in list_folders:
-                if os.path.isdir(folder):
-                    shutil.rmtree(folder)
-            for extension in list_file_extensions:
-                os.remove(os.path.join(file_path, case + extension))
-        pass
+        folders = ['hale/output']
+        for folder in folders:
+            shutil.rmtree(self.route_file_dir + '/' + folder)
+        files = ['hale/hale.aero.h5', 'hale/hale.fem.h5', 'hale/hale.sharpy']
+        for file in files:
+            file_dir = self.route_file_dir + '/' + file
+            if os.path.isfile(file_dir):
+                os.remove(file_dir)
+            
+   
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file
diff --git a/tests/sourcepanelmethod/__init__.py b/tests/sourcepanelmethod/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/sourcepanelmethod/geometry_parameter_models.json b/tests/sourcepanelmethod/geometry_parameter_models.json
new file mode 100644
index 000000000..06aae741b
--- /dev/null
+++ b/tests/sourcepanelmethod/geometry_parameter_models.json
@@ -0,0 +1 @@
+{"phantom_wing": {"length_radius_ratio": 50.0, "radius_chord_ratio": 0.2, "vertical_wing_position": 0.0, "radius_half_wingspan_ratio": 0.02, "length_offset_nose_to_wing_ratio": 0.5, "fuselage_shape": "cylindrical"}, "low_wing": {"length_radius_ratio": 16.0, "radius_chord_ratio": 0.5, "vertical_wing_position": -0.5, "radius_half_wingspan_ratio": 0.222, "fuselage_shape": "cylindrical", "length_offset_nose_to_wing_ratio": 0.5}}
\ No newline at end of file
diff --git a/tests/sourcepanelmethod/test_data/results_low_wing_coupled_0.csv b/tests/sourcepanelmethod/test_data/results_low_wing_coupled_0.csv
new file mode 100644
index 000000000..d7cd3e6e9
--- /dev/null
+++ b/tests/sourcepanelmethod/test_data/results_low_wing_coupled_0.csv
@@ -0,0 +1,21 @@
+0.000000000000000000e+00 0.000000000000000000e+00
+0.000000000000000000e+00 0.000000000000000000e+00
+0.000000000000000000e+00 0.000000000000000000e+00
+0.000000000000000000e+00 0.000000000000000000e+00
+5.630631000000000386e-01 2.637046000000000112e-01
+7.038288000000000322e-01 2.587818000000000063e-01
+8.445945999999999732e-01 2.560145999999999811e-01
+9.853604000000000251e-01 2.521720000000000073e-01
+1.126125999999999960e+00 2.474606000000000028e-01
+1.266891999999999907e+00 2.419823000000000113e-01
+1.407658000000000076e+00 2.357548999999999895e-01
+1.548423000000000105e+00 2.287254999999999983e-01
+1.689189000000000052e+00 2.207763000000000086e-01
+1.829954999999999998e+00 2.117207000000000117e-01
+1.970720999999999945e+00 2.012865000000000071e-01
+2.111486000000000196e+00 1.890803000000000067e-01
+2.252251999999999921e+00 1.745153000000000121e-01
+2.393018000000000089e+00 1.566525000000000001e-01
+2.533783999999999814e+00 1.337860999999999911e-01
+2.674549999999999983e+00 1.019361000000000017e-01
+2.815315000000000012e+00 8.357517000000000418e-02
diff --git a/tests/sourcepanelmethod/test_data/results_low_wing_coupled_1.csv b/tests/sourcepanelmethod/test_data/results_low_wing_coupled_1.csv
new file mode 100644
index 000000000..f558b3649
--- /dev/null
+++ b/tests/sourcepanelmethod/test_data/results_low_wing_coupled_1.csv
@@ -0,0 +1,21 @@
+0.000000000000000000e+00 0.000000000000000000e+00
+0.000000000000000000e+00 0.000000000000000000e+00
+0.000000000000000000e+00 0.000000000000000000e+00
+0.000000000000000000e+00 0.000000000000000000e+00
+5.630631000000000386e-01 2.636864999999999903e-01
+7.038288000000000322e-01 2.592814999999999981e-01
+8.445945999999999732e-01 2.565006000000000230e-01
+9.853604000000000251e-01 2.526663000000000103e-01
+1.126125999999999960e+00 2.479610999999999899e-01
+1.266891999999999907e+00 2.424882999999999900e-01
+1.407658000000000076e+00 2.362653000000000114e-01
+1.548423000000000105e+00 2.292386999999999897e-01
+1.689189000000000052e+00 2.212899000000000116e-01
+1.829954999999999998e+00 2.122312000000000087e-01
+1.970720999999999945e+00 2.017887999999999904e-01
+2.111486000000000196e+00 1.895677000000000056e-01
+2.252251999999999921e+00 1.749785000000000090e-01
+2.393018000000000089e+00 1.570786999999999878e-01
+2.533783999999999814e+00 1.341562000000000032e-01
+2.674549999999999983e+00 1.022175000000000028e-01
+2.815315000000000012e+00 8.376727999999999952e-02
diff --git a/tests/sourcepanelmethod/test_source_panel_method.py b/tests/sourcepanelmethod/test_source_panel_method.py
new file mode 100644
index 000000000..7056e6475
--- /dev/null
+++ b/tests/sourcepanelmethod/test_source_panel_method.py
@@ -0,0 +1,115 @@
+import unittest
+import os
+import numpy as np
+from sharpy.cases.templates.fuselage_wing_configuration.fuselage_wing_configuration import Fuselage_Wing_Configuration
+from sharpy.cases.templates.fuselage_wing_configuration.fwc_get_settings import define_simulation_settings
+
+
+
+class TestSourcePanelMethod(unittest.TestCase):
+
+    route_test_dir = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))
+
+    def test_ellipsoid(self):
+        """
+        Computes the pressure distribution over an ellipsoid. The results should
+        match the analyitcal solution according to potential flow theory (see Chapter 5
+        Katz, J., and Plotkin, A., Low-speed aerodynamics, Vol. 13, Cambridge University 
+        Press, 2001).
+        """
+
+        # define model variables
+        radius_ellipsoid = 0.2
+        length_ellipsoid = 2.
+        u_inf = 10 
+        alpha_deg = 0
+        n_elem = 30
+        num_radial_panels  = 24
+        lifting_only= False
+        fuselage_shape = 'ellipsoid'
+        fuselage_discretisation = 'uniform'
+        # define case name and folders
+        case_route = self.route_test_dir + '/cases/'
+        output_route = self.route_test_dir + '/output/'
+        case_name = 'ellipsoid'
+        enforce_uniform_fuselage_discretisation = True
+
+        # generate ellipsoid model
+        ellipsoidal_body = Fuselage_Wing_Configuration(case_name, case_route, output_route) 
+        ellipsoidal_body.init_aeroelastic(lifting_only=lifting_only,
+                                          max_radius=radius_ellipsoid,
+                                      fuselage_length=length_ellipsoid,
+                                      offset_nose_wing=length_ellipsoid/2,
+                                      n_elem_fuselage=n_elem,
+                                      num_radial_panels=num_radial_panels,
+                                      fuselage_shape=fuselage_shape,
+                                      enforce_uniform_fuselage_discretisation=enforce_uniform_fuselage_discretisation,
+                                      fuselage_discretisation=fuselage_discretisation)
+        ellipsoidal_body.generate()
+        
+        # define settings
+        flow = ['BeamLoader',
+                'NonliftingbodygridLoader',
+                'StaticUvlm',
+                'WriteVariablesTime'
+                ]
+        settings = define_simulation_settings(flow, 
+                                              ellipsoidal_body, 
+                                              alpha_deg, 
+                                              u_inf, 
+                                              lifting_only=False, 
+                                              nonlifting_only=True, 
+                                              horseshoe=True,
+                                              writeCpVariables=True)
+        ellipsoidal_body.create_settings(settings)
+
+        # run simulation
+        ellipsoidal_body.run()
+
+        # postprocess
+        cp_distribution_SHARPy = self.load_pressure_distribution(output_route + '/' + case_name + '/WriteVariablesTime/', 
+                                                                 ellipsoidal_body.structure.n_node_fuselage)
+        dx = length_ellipsoid/(ellipsoidal_body.structure.n_node_fuselage-1)
+        x_collocation_points = np.linspace(-length_ellipsoid/2+dx/2, length_ellipsoid/2-dx/2, ellipsoidal_body.structure.n_node_fuselage)
+        cp_distribution_analytcal = self.get_analytical_pressure_distribution(radius_ellipsoid, x_collocation_points)
+
+        # check results
+        with self.subTest('pressure_coefficient'):
+            # Higher deviations near the nose and tail due to local coarser discretisation. Check only mid-body values!
+            np.testing.assert_array_almost_equal(cp_distribution_SHARPy[4:-4], cp_distribution_analytcal[4:-4], decimal=3)
+
+    def load_pressure_distribution(self, output_folder, n_collocation_points):
+        """
+            Loads the resulting pressure coefficients saved in txt-files.
+        """
+        cp_distribution = np.zeros((n_collocation_points,))
+        for i_collocation_point in range(n_collocation_points):
+            cp_distribution[i_collocation_point] = np.loadtxt(output_folder + 'nonlifting_pressure_coefficients_panel_isurf0_im0_in{}.dat'.format(i_collocation_point))[1]
+        return cp_distribution
+
+    def get_analytical_pressure_distribution(self, radius, x_coordinates):
+        """
+            Computes the analytical solution of the pressure distribution over
+            an ellipsoid in potential flow for the previous specified ellipsoid
+            model. Equations used are taken from
+                https://www.symscape.com/examples/panel/potential_flow_ellipsoid
+        """
+        a = np.sqrt(1 - radius**2)
+        b = 2 * ((1-a**2)/a**3) * (np.arctanh(a)-a)
+        u = 2./(2.-b) * np.sqrt((1-x_coordinates**2)/(1-x_coordinates**2 * a**2))
+        return 1-u**2
+    
+    def tearDown(self):
+        """
+            Removes all created files within this test.
+        """
+        import shutil
+        folders = ['cases', 'output']
+        for folder in folders:
+            shutil.rmtree(self.route_test_dir + '/' + folder)
+
+
+if __name__ == '__main__':
+    import unittest
+
+    unittest.main()
diff --git a/tests/sourcepanelmethod/test_vlm_coupled_spm.py b/tests/sourcepanelmethod/test_vlm_coupled_spm.py
new file mode 100644
index 000000000..b5b651799
--- /dev/null
+++ b/tests/sourcepanelmethod/test_vlm_coupled_spm.py
@@ -0,0 +1,296 @@
+import unittest
+import os
+import numpy as np
+from sharpy.cases.templates.fuselage_wing_configuration.fuselage_wing_configuration import Fuselage_Wing_Configuration
+from sharpy.cases.templates.fuselage_wing_configuration.fwc_get_settings import define_simulation_settings
+import json
+
+class TestUvlmCoupledWithSourcePanelMethod(unittest.TestCase):
+
+ 
+    def test_phantom_panels(self):
+        """
+            Phantom Panel Test (Steady and Dynamic)
+            The lift distribution over a rectangular high-aspect ratio
+            wing is computed. First a wing_only configuration is considered,
+            while second, we activate the phantom panels created within the 
+            fuselage although, the effect of the source panels is omitted. 
+            With the defined interpolation scheme, the same lift distribution 
+            must be obtained.
+        """
+        self.define_folder()
+        
+
+        model = 'phantom_wing'
+        fuselage_length = 10
+        dict_geometry_parameters = self.get_geometry_parameters(model,
+                                                                fuselage_length=fuselage_length)
+
+        # Freestream Conditions
+        alpha_deg = 5
+        u_inf = 10
+
+        # Discretization
+        dict_discretization = {
+            'n_elem_per_wing': 20,
+            'n_elem_fuselage': 10,
+            'num_chordwise_panels': 4
+        }
+
+        # Simulation settings
+        horseshoe = True
+        list_phantom_test = [False, True]
+        list_dynamic_test = [False, True]
+        dynamic_structural_solver = 'NonLinearDynamicPrescribedStep'
+
+        # Numerical Settings for DynamicCoupled (very small for faster test runs)
+        fsi_substeps = 2
+        fsi_tolerance = 1e-2
+
+        # Simlation Solver Flow
+        flow = ['BeamLoader',
+                'AerogridLoader',
+                'NonliftingbodygridLoader',
+                'StaticUvlm',
+                'StaticCoupled',
+                'BeamLoads',
+                'LiftDistribution',
+                'DynamicCoupled',
+                'LiftDistribution',
+                    ]
+        # define model variables
+        for dynamic in list_dynamic_test:
+            list_results_lift_distribution = []
+            for phantom_test in list_phantom_test:
+                horseshoe = not dynamic
+                lifting_only = not phantom_test
+                case_name = model + '_phantom{}_dynamic{}'.format(int(phantom_test),
+                                                                  int(dynamic))
+                
+                # generate ellipsoid model
+                phantom_wing = self.generate_model(case_name, 
+                                                dict_geometry_parameters,
+                                                dict_discretization, 
+                                                lifting_only)
+                # Adjust flow for case
+                flow_case = flow.copy()
+
+                if lifting_only:
+                    n_tsteps = 5
+                    flow_case.remove('NonliftingbodygridLoader')
+                if not dynamic:
+                    n_tsteps = 0
+                    flow_case.remove('StaticCoupled')
+                    flow_case.remove('DynamicCoupled')
+                self.generate_simulation_settings(flow_case, 
+                                                phantom_wing, 
+                                                alpha_deg, 
+                                                u_inf, 
+                                                lifting_only,
+                                                n_tsteps = n_tsteps,
+                                                horseshoe=horseshoe,
+                                                phantom_test=phantom_test,
+                                                dynamic_structural_solver=dynamic_structural_solver,
+                                                fsi_substeps=fsi_substeps,
+                                                fsi_tolerance=fsi_tolerance)
+                # # run simulation
+                phantom_wing.run()
+
+                # get results
+                list_results_lift_distribution.append(self.load_lift_distribution(
+                    self.output_route + '/' + case_name,
+                    phantom_wing.structure.n_node_right_wing,
+                    n_tsteps
+                )) 
+            
+            # check results
+            with self.subTest(self.get_test_name('phantom', dynamic)):
+                np.testing.assert_array_almost_equal(list_results_lift_distribution[0][3:, 1], 
+                                                     list_results_lift_distribution[1][3:, 1], 
+                                                     decimal=2) #3 - int(dynamic))
+
+
+    
+    def test_fuselage_wing_configuration(self):
+        """
+            Test spanwise lift distribution on a fuselage-wing configuration.
+            
+            The lift distribution on low wing configuration is computed considering
+            fuselage effects. The final results are compared to a previous solution 
+            (backward compatibility) that matches the experimental lift distribution for this case. 
+        """
+        self.define_folder()
+        model = 'low_wing'
+        fuselage_length = 10
+        dict_geometry_parameters = self.get_geometry_parameters(model,
+                                                                fuselage_length=fuselage_length)
+
+        # Freestream Conditions
+        alpha_deg = 2.9
+        u_inf = 10
+
+        # Discretization
+        dict_discretization = {
+            'n_elem_per_wing': 10,
+            'n_elem_fuselage': 30,
+            'num_chordwise_panels': 8,
+            'num_radial_panels': 36
+        }
+
+        # Simulation settings
+        horseshoe = True
+        phantom_test = False
+        lifting_only = False
+        # Simlation Solver Flow
+        flow = ['BeamLoader',
+                'AerogridLoader',
+                'NonliftingbodygridLoader',
+                'StaticUvlm',
+                'StaticCoupled',
+                'BeamLoads',
+                'LiftDistribution',
+                    ]
+        for static_coupled_solver in [False, True]:
+            case_name = '{}_coupled_{}'.format(model, int(static_coupled_solver))
+            wing_fuselage_model = self.generate_model(case_name, 
+                                                dict_geometry_parameters,
+                                                dict_discretization, 
+                                                lifting_only)
+
+            flow_case = flow.copy()
+            if static_coupled_solver:
+                flow_case.remove('StaticUvlm')
+            else:
+                flow_case.remove('StaticCoupled')
+            self.generate_simulation_settings(flow_case, 
+                                                wing_fuselage_model, 
+                                                alpha_deg, 
+                                                u_inf, 
+                                                lifting_only,
+                                                horseshoe=horseshoe,
+                                                phantom_test=phantom_test)
+            # run simulation
+            wing_fuselage_model.run()
+            # get results
+            lift_distribution = self.load_lift_distribution(
+                self.output_route + case_name,
+                wing_fuselage_model.structure.n_node_right_wing
+                )
+            # check results
+            lift_distribution_test = np.loadtxt(self.route_test_dir + "/test_data/results_{}.csv".format(case_name))
+            with self.subTest(self.get_test_name('lift distribution and spanwise wing deformation',
+                                                 static_coupled=static_coupled_solver)):        
+                np.testing.assert_array_almost_equal(lift_distribution_test, lift_distribution, decimal=3)
+
+    def get_test_name(self, base_name, dynamic=False, static_coupled=False):
+        if static_coupled:
+            base_name += ' coupled'
+        else:
+            if dynamic:
+                base_name += ' uvlm'
+            else:
+                base_name += ' vlm'
+        return base_name
+    def get_geometry_parameters(self, model_name,fuselage_length=10):
+        """
+            Geometry parameters are loaded from json init file for the specified model. 
+            Next, final geoemtry parameteres, depending on the fuselage length are 
+            calculated and return within a dict.
+        """
+        with open(self.route_test_dir + '/geometry_parameter_models.json', 'r') as fp:
+            parameter_models = json.load(fp)[model_name]
+
+        geometry_parameters = {
+            'fuselage_length': fuselage_length,         
+            'max_radius': fuselage_length/parameter_models['length_radius_ratio'],
+            'fuselage_shape': parameter_models['fuselage_shape'],
+        }
+        geometry_parameters['chord']=geometry_parameters['max_radius']/parameter_models['radius_chord_ratio']
+        geometry_parameters['half_wingspan'] = geometry_parameters['max_radius']/parameter_models['radius_half_wingspan_ratio']
+        geometry_parameters['offset_nose_wing'] = parameter_models['length_offset_nose_to_wing_ratio'] * fuselage_length
+        geometry_parameters['vertical_wing_position'] = parameter_models['vertical_wing_position'] * geometry_parameters['max_radius']
+        return geometry_parameters
+
+    def generate_model(self, 
+                       case_name, 
+                       dict_geometry_parameters,
+                       dict_discretisation, 
+                       lifting_only):
+        """
+            Aircraft model object is generated and structural and aerodynamic (lifting and nonlifting)
+            input files are generated.
+        """
+        aircraft_model = Fuselage_Wing_Configuration(case_name, self.case_route, self.output_route)
+        aircraft_model.init_aeroelastic(lifting_only=lifting_only,
+                                    **dict_discretisation,
+                                    **dict_geometry_parameters)
+        aircraft_model.generate()
+        return aircraft_model
+    
+    def generate_simulation_settings(self, 
+                                     flow, 
+                                     aircraft_model, 
+                                     alpha_deg, 
+                                     u_inf, 
+                                     lifting_only,
+                                     n_tsteps=0,
+                                     horseshoe=True,
+                                     nonlifting_only=False,
+                                     phantom_test=False,
+                                     dynamic_structural_solver='NonLinearDynamicPrescribedStep',
+                                     fsi_substeps=200,
+                                     fsi_tolerance=1e-6
+                                     ):
+        """
+            Simulation settings are defined and written to the ".sharpy" input file.
+        """
+        settings = define_simulation_settings(flow, 
+                                                aircraft_model, 
+                                                alpha_deg, 
+                                                u_inf, 
+                                                dt=self.get_timestep(aircraft_model, u_inf),
+                                                n_tsteps=n_tsteps,
+                                                lifting_only=lifting_only, 
+                                                phantom_test=phantom_test, 
+                                                nonlifting_only=nonlifting_only, 
+                                                horseshoe=horseshoe,
+                                                dynamic_structural_solver=dynamic_structural_solver,
+                                                fsi_substeps=fsi_substeps,
+                                                fsi_tolerance=fsi_tolerance)
+        aircraft_model.create_settings(settings)
+
+    def get_timestep(self, model, u_inf):
+        return model.aero.chord_wing / model.aero.num_chordwise_panels / u_inf
+    
+    def define_folder(self):
+        """
+            Initializes all folder path needed and creates case folder.
+        """
+        self.route_test_dir = os.path.abspath(os.path.dirname(os.path.realpath(__file__)))
+        self.case_route = self.route_test_dir + '/cases/'
+        self.output_route = self.route_test_dir + '/output/'
+    
+        if not os.path.exists(self.case_route):
+            os.makedirs(self.case_route)
+        
+    def load_lift_distribution(self, output_folder, n_node_wing, n_tsteps=0):
+        """
+            Loads the resulting pressure coefficients saved in txt-files.
+        """
+        lift_distribution = np.loadtxt(output_folder + '/liftdistribution/liftdistribution_ts{}.txt'.format(str(n_tsteps)), delimiter=',')
+        return lift_distribution[:n_node_wing,[1,-1]]
+
+    def tearDown(self):
+        """
+            Removes all created files within this test.
+        """
+        import shutil
+        folders = ['cases', 'output']
+        for folder in folders:
+            shutil.rmtree(self.route_test_dir + '/' + folder)
+
+
+if __name__ == '__main__':
+    import unittest
+
+    unittest.main()