diff --git a/sharpy/generators/polaraeroforces.py b/sharpy/generators/polaraeroforces.py
index 9670690e5..8f2a5c6f2 100644
--- a/sharpy/generators/polaraeroforces.py
+++ b/sharpy/generators/polaraeroforces.py
@@ -9,7 +9,7 @@
 
 @generator_interface.generator
 class PolarCorrection(generator_interface.BaseGenerator):
-    """
+    r"""
     This generator corrects the aerodynamic forces from UVLM based on the airfoil polars provided by the user in the
     ``aero.h5`` file. Polars are entered for each airfoil, in a table comprising ``AoA (rad), CL, CD, CM``.
 
diff --git a/sharpy/linear/assembler/linearbeam.py b/sharpy/linear/assembler/linearbeam.py
index 8e9c8e50c..311906a68 100644
--- a/sharpy/linear/assembler/linearbeam.py
+++ b/sharpy/linear/assembler/linearbeam.py
@@ -370,7 +370,7 @@ def remove_symmetric_modes(self):
         self.sys.num_modes = len(self.sys.freq_natural)
 
     def unpack_ss_vector(self, x_n, u_n, struct_tstep):
-        """
+        r"""
         Warnings:
             Under development. Missing:
                 * Accelerations
diff --git a/sharpy/linear/assembler/lineargustassembler.py b/sharpy/linear/assembler/lineargustassembler.py
index 72b6878ca..c03e0ef73 100644
--- a/sharpy/linear/assembler/lineargustassembler.py
+++ b/sharpy/linear/assembler/lineargustassembler.py
@@ -205,7 +205,7 @@ def assemble_gust_statespace(self, a_i, b_i, c_i, d_i):
         return gustss
 
     def discretise_domain(self):
-        """
+        r"""
         Generates a "gust-station" domain, aligned with the free-stream velocity equispaced in
         :math:`\Delta t U_\infty` increments.
 
diff --git a/sharpy/postproc/aerogridplot.py b/sharpy/postproc/aerogridplot.py
index 88a7f2b97..b78b56728 100644
--- a/sharpy/postproc/aerogridplot.py
+++ b/sharpy/postproc/aerogridplot.py
@@ -1,7 +1,6 @@
 import os
 
 import numpy as np
-from tvtk.api import tvtk, write_data
 
 import sharpy.utils.algebra as algebra
 import sharpy.utils.cout_utils as cout
@@ -11,6 +10,8 @@
 import sharpy.aero.utils.uvlmlib as uvlmlib
 from sharpy.utils.constants import vortex_radius_def
 
+from sharpy.utils.plotutils import plot_frame_to_vtk
+
 
 @solver
 class AerogridPlot(BaseSolver):
@@ -131,7 +132,7 @@ def initialise(self, data, custom_settings=None, caller=None, restart=False):
 
     def run(self, **kwargs):
     
-        online = su.set_value_or_default(kwargs, 'online', False)
+        online = kwargs.get('online', False)
 
         # TODO: Create a dictionary to plot any variable as in beamplot
         if not online:
@@ -167,205 +168,88 @@ def plot_body(self):
                                                    struct_tstep)
 
         for i_surf in range(aero_tstep.n_surf):
-            filename = (self.body_filename +
-                        '_' +
-                        ('%02u_' % i_surf) +
-                        ('%06u' % self.ts) +
-                        '.vtu')
+            filename = f"{self.body_filename}_{i_surf:02d}_{self.ts:06d}.vtu"
+
+            zeta = np.moveaxis(aero_tstep.zeta[i_surf], 0, -1).copy()  # [m+1, n+1, 3]
+            zeta_dot = np.moveaxis(aero_tstep.zeta_dot[i_surf][:3, ...], 0, -1).copy()  # [3, m+1, n+1]
+            gamma = aero_tstep.gamma[i_surf].copy()  # [m, n]
+            gamma_dot = aero_tstep.gamma_dot[i_surf].copy()  # [m, n]
+            normals = np.moveaxis(aero_tstep.normals[i_surf], 0, -1).copy()  # [m, n, 3]
+            surf_id = np.ones(gamma.shape, dtype=int) * i_surf
+
+            if self.settings['include_rbm']:
+                if struct_tstep.mb_dict is None:
+                    zeta += np.expand_dims(struct_tstep.for_pos[:3], (0, 1))
+                else:
+                    zeta += np.expand_dims(self.data.structure.timestep_info[0].cga() @ struct_tstep.for_pos[:3], (0, 1))
+            if self.settings['include_forward_motion']:
+                zeta[..., 0] -= self.settings['dt'] * self.ts * self.settings['u_inf']
 
-            dims = aero_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_gamma = np.zeros((panel_data_dim,))
-            panel_gamma_dot = 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))
-            point_unsteady_cf = np.zeros((point_data_dim, 3))
-            zeta_dot = np.zeros((point_data_dim, 3))
-            u_inf = np.zeros((point_data_dim, 3))
-            if self.settings['include_velocities']:
-                vel = np.zeros((point_data_dim, 3))
-            counter = -1
+            if self.settings['include_incidence_angle']:
+                incidence_angle = aero_tstep.postproc_cell['incidence_angle'][i_surf]
+            else:
+                incidence_angle = np.zeros(gamma.shape)
 
-            # coordinates of corners
-            for i_n in range(dims[1]+1):
-                for i_m in range(dims[0]+1):
-                    counter += 1
-                    coords[counter, :] = aero_tstep.zeta[i_surf][:, i_m, i_n]
-                    if self.settings['include_rbm']:
-                        #TODO: fix for lack of g frame description in nonlineardynamicmultibody.py
-                        if struct_tstep.mb_dict is None:
-                            coords[counter, :] += struct_tstep.for_pos[0:3]
-                        else:
-                            #TODO: uncomment for dynamic trim
-                            # try:
-                            #     cga = algebra.euler2rot([0, self.data.trimmed_values[0], 0])
-                            #     coords[counter, :] += np.dot(cga, struct_tstep.for_pos[0:3])
-                            # except AttributeError:
-                            coords[counter, :] += np.dot(self.data.structure.timestep_info[0].cga(), struct_tstep.for_pos[0:3])
-                    if self.settings['include_forward_motion']:
-                        coords[counter, 0] -= self.settings['dt']*self.ts*self.settings['u_inf']
 
-            if self.settings['include_incidence_angle']:
-                incidence_angle = np.zeros_like(panel_gamma)
+            point_struct_id = np.broadcast_to(np.expand_dims(self.data.aero.aero2struct_mapping[i_surf], 0), zeta.shape[:2])
+            point_cf = np.moveaxis(aero_tstep.forces[i_surf][:3, ...], 0, -1)
 
-            counter = -1
-            node_counter = -1
-            for i_n in range(dims[1] + 1):
-                global_counter = self.data.aero.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, :] = aero_tstep.forces[i_surf][0:3, i_m, i_n]
-                    try:
-                        point_unsteady_cf[node_counter, :] = aero_tstep.dynamic_forces[i_surf][0:3, i_m, i_n]
-                    except AttributeError:
-                        pass
-                    try:
-                        zeta_dot[node_counter, :] = aero_tstep.zeta_dot[i_surf][0:3, i_m, i_n]
-                    except AttributeError:
-                        pass
-                    try:
-                        u_inf[node_counter, :] = aero_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
+            try:
+                point_unsteady_cf = np.moveaxis(aero_tstep.dynamic_forces[i_surf][:3, ...], 0, -1)
+            except AttributeError:
+                point_unsteady_cf = np.zeros(zeta.shape)
 
-                    conn.append([node_counter + 0,
-                                 node_counter + 1,
-                                 node_counter + dims[0]+2,
-                                 node_counter + dims[0]+1])
-                    # cell data
-                    normal[counter, :] = aero_tstep.normals[i_surf][:, i_m, i_n]
-                    panel_id[counter] = counter
-                    panel_surf_id[counter] = i_surf
-                    panel_gamma[counter] = aero_tstep.gamma[i_surf][i_m, i_n]
-                    panel_gamma_dot[counter] = aero_tstep.gamma_dot[i_surf][i_m, i_n]
+            try:
+                u_inf = np.moveaxis(aero_tstep.u_ext[i_surf][:3, ...], 0, -1)
+            except AttributeError:
+                u_inf = np.zeros(zeta.shape)
 
-                    if self.settings['include_incidence_angle']:
-                        incidence_angle[counter] = \
-                            aero_tstep.postproc_cell['incidence_angle'][i_surf][i_m, i_n]
 
             if self.settings['include_velocities']:
                 vel = uvlmlib.uvlm_calculate_total_induced_velocity_at_points(aero_tstep,
-                                                                              coords,
+                                                                              zeta.reshape(-1, 3),
                                                                               self.settings['vortex_radius'],
                                                                               struct_tstep.for_pos,
-                                                                              self.settings['num_cores'])
+                                                                              self.settings['num_cores']).reshape(zeta.shape)
+            else:
+                vel = np.zeros(zeta.shape)
+
+            plot_frame_to_vtk(zeta,
+                              filename,
+                              node_scalar_data={'point_struct_id': point_struct_id},
+                              node_vector_data={'zeta_dot': zeta_dot, 'point_steady_force': point_cf, 'point_unsteady_force': point_unsteady_cf, 'u_inf': u_inf, 'velocity': vel},
+                              cell_scalar_data={'gamma': gamma, "gamma_dot": gamma_dot, "incidence_angle": incidence_angle, 'surf_id': surf_id},
+                              cell_vector_data={'panel_normal': normals},)
 
-            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_gamma)
-            ug.cell_data.get_array(2).name = 'panel_gamma'
-            ug.cell_data.add_array(panel_gamma_dot)
-            ug.cell_data.get_array(3).name = 'panel_gamma_dot'
-            if self.settings['include_incidence_angle']:
-                ug.cell_data.add_array(incidence_angle)
-                ug.cell_data.get_array(4).name = 'incidence_angle'
-            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(point_unsteady_cf)
-            ug.point_data.get_array(3).name = 'point_unsteady_force'
-            ug.point_data.add_array(zeta_dot)
-            ug.point_data.get_array(4).name = 'zeta_dot'
-            ug.point_data.add_array(u_inf)
-            ug.point_data.get_array(5).name = 'u_inf'
-            if self.settings['include_velocities']:
-                ug.point_data.add_array(vel)
-                ug.point_data.get_array(6).name = 'velocity'
-            write_data(ug, filename)
 
     def plot_wake(self):
         for i_surf in range(self.data.aero.timestep_info[self.ts].n_surf):
-            filename = (self.wake_filename +
-                        '_' +
-                        ('%02u_' % i_surf) +
-                        ('%06u' % self.ts) +
-                        '.vtu')
-
-            dims_star = self.data.aero.timestep_info[self.ts].dimensions_star[i_surf, :].copy()
-            dims_star[0] -= self.settings['minus_m_star']
+            filename = f"{self.wake_filename}_{i_surf:02d}_{self.ts:06d}.vtu"
 
-            point_data_dim = (dims_star[0]+1)*(dims_star[1]+1)
-            panel_data_dim = (dims_star[0])*(dims_star[1])
+            zeta_star = np.moveaxis(self.data.aero.timestep_info[self.ts].zeta_star[i_surf], 0, -1).copy() # [m_star+1, n+1, 3]
+            gamma = self.data.aero.timestep_info[self.ts].gamma_star[i_surf].copy()  # [m_star, n]
 
-            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_gamma = np.zeros((panel_data_dim,))
-            counter = -1
-            # rotation_mat = self.data.structure.timestep_info[self.ts].cga().T
-            # coordinates of corners
-            for i_n in range(dims_star[1]+1):
-                for i_m in range(dims_star[0]+1):
-                    counter += 1
-                    coords[counter, :] = self.data.aero.timestep_info[self.ts].zeta_star[i_surf][:, i_m, i_n]
-                    if self.settings['include_rbm']:
-                        #TODO: fix for lack of g frame description in nonlineardynamicmultibody.py
-                        if self.data.structure.timestep_info[self.ts].mb_dict is None:
-                            coords[counter, :] += self.data.structure.timestep_info[self.ts].for_pos[0:3]
-                        else:
-                            #TODO: uncomment for dynamic trim
-                            # try:
-                            #     cga = algebra.euler2rot([0, self.data.trimmed_values[0], 0])
-                            #     coords[counter, :] += np.dot(cga, self.data.structure.timestep_info[self.ts].for_pos[0:3])
-                            # except AttributeError:
-                            coords[counter, :] += np.dot(self.data.structure.timestep_info[0].cga(), self.data.structure.timestep_info[self.ts].for_pos[0:3])                                                        
-                    if self.settings['include_forward_motion']:
-                        coords[counter, 0] -= self.settings['dt']*self.ts*self.settings['u_inf']
+            surf_id = np.ones(gamma.shape, dtype=int) * i_surf
 
-            counter = -1
-            node_counter = -1
-            # wake
-            for i_n in range(dims_star[1]+1):
-                for i_m in range(dims_star[0]+1):
-                    node_counter += 1
-                    # cell data
-                    if i_n < dims_star[1] and i_m < dims_star[0]:
-                        counter += 1
-                    else:
-                        continue
+            if self.settings['include_rbm']:
+                if self.data.structure.timestep_info[self.ts].mb_dict is None:
+                    zeta_star += np.expand_dims(self.data.structure.timestep_info[self.ts].for_pos[:3], (0, 1))
+                else:
+                    zeta_star += np.expand_dims(self.data.structure.timestep_info[0].cga()
+                                                @ self.data.structure.timestep_info[self.ts].for_pos[:3], (0, 1))
+            if self.settings['include_forward_motion']:
+                zeta_star[..., 0] -= self.settings['dt'] * self.ts * self.settings['u_inf']
 
-                    conn.append([node_counter + 0,
-                                 node_counter + 1,
-                                 node_counter + dims_star[0]+2,
-                                 node_counter + dims_star[0]+1])
-                    panel_id[counter] = counter
-                    panel_surf_id[counter] = i_surf
-                    panel_gamma[counter] = self.data.aero.timestep_info[self.ts].gamma_star[i_surf][i_m, i_n]
+            plot_frame_to_vtk(
+                zeta_star,
+                filename,
+                node_vector_data={},
+                cell_scalar_data={
+                    "gamma": gamma,
+                    "surf_id": surf_id,
+                },
+            )
 
-            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_gamma)
-            ug.cell_data.get_array(2).name = 'panel_gamma'
-            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]
@@ -450,21 +334,3 @@ def plot_nonlifting_surfaces(self):
             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/beamplot.py b/sharpy/postproc/beamplot.py
index 0e21c0763..6c36c57f3 100644
--- a/sharpy/postproc/beamplot.py
+++ b/sharpy/postproc/beamplot.py
@@ -1,7 +1,6 @@
 import os
 
 import numpy as np
-from tvtk.api import tvtk, write_data
 
 import sharpy.utils.cout_utils as cout
 from sharpy.utils.solver_interface import solver, BaseSolver
@@ -9,6 +8,11 @@
 import sharpy.utils.algebra as algebra
 
 
+import vtk
+from vtk.numpy_interface import algorithms as algs
+from vtk.numpy_interface import dataset_adapter as dsa
+
+
 @solver
 class BeamPlot(BaseSolver):
     """
@@ -118,9 +122,7 @@ def plot(self, online):
                 self.write_for(it)
 
     def write_beam(self, it):
-        it_filename = (self.filename +
-                       ('%06u' % it) +
-                       '.vtu')
+        it_filename = f"{self.filename}{it:06d}.vtp"
         num_nodes = self.data.structure.num_node
         num_elem = self.data.structure.num_elem
 
@@ -132,7 +134,6 @@ def write_beam(self, it):
         local_y = np.zeros((num_nodes, 3))
         local_z = np.zeros((num_nodes, 3))
         coords_a = np.zeros((num_nodes, 3))
-
         app_forces = np.zeros((num_nodes, 3))
         app_moment = np.zeros((num_nodes, 3))
 
@@ -150,14 +151,9 @@ def write_beam(self, it):
         if tstep.mb_dict is None:
             pass
         else:
-            #TODO: fix for lack of g frame description in nonlineardynamicmultibody.py
             for i_node in range(tstep.num_node):
-                #TODO: uncomment for dynamic trim
-                # try:
-                #     c = algebra.euler2rot([0, self.data.trimmed_values[0], 0])
-                # except AttributeError:
                 c = self.data.structure.timestep_info[0].cga()
-                coords[i_node, :] += np.dot(c, tstep.for_pos[0:3])
+                coords[i_node, :] += c @ tstep.for_pos[:3]
 
         # check if I can output gravity forces
         with_gravity = False
@@ -260,115 +256,128 @@ def write_beam(self, it):
             conn[i_elem, :] = self.data.structure.elements[i_elem].reordered_global_connectivities
             elem_id[i_elem] = i_elem
 
-        ug = tvtk.UnstructuredGrid(points=coords)
-        ug.set_cells(tvtk.Line().cell_type, conn)
-        ug.cell_data.scalars = elem_id
-        ug.cell_data.scalars.name = 'elem_id'
-        counter = 1
+        ug = vtk.vtkPolyData(points=coords)
+        cells = vtk.vtkCellArray()
+        for _conn in conn:
+            line = vtk.vtkPolyLine()
+            line.GetPointIds().SetNumberOfIds(3)
+            line.GetPointIds().SetId(0, _conn[0])
+            line.GetPointIds().SetId(1, _conn[1])
+            line.GetPointIds().SetId(2, _conn[2])
+            cells.InsertNextCell(line)
+        ug.SetLines(cells)
+
+        ug.GetCellData().AddArray(dsa.numpyTovtkDataArray(np.array(elem_id), name='elem_id'))
+        ug.GetCellData().AddArray(dsa.numpyTovtkDataArray(coords_a_cell, name="coords_a_elem"))
+        ug.GetPointData().AddArray(
+            dsa.numpyTovtkDataArray(node_id, name="node_id")
+        )
+        ug.GetPointData().AddArray(dsa.numpyTovtkDataArray(local_x, name="local_x"))
+        ug.GetPointData().AddArray(dsa.numpyTovtkDataArray(local_y, name="local_y"))
+        ug.GetPointData().AddArray(dsa.numpyTovtkDataArray(local_z, name="local_z"))
+        ug.GetPointData().AddArray(dsa.numpyTovtkDataArray(coords_a, name="coords_a"))
+
         if with_postproc_cell:
             for k in postproc_cell_vector:
-                ug.cell_data.add_array(tstep.postproc_cell[k])
-                ug.cell_data.get_array(counter).name = k + '_cell'
-                counter += 1
+                ug.GetCellData().AddArray(
+                    dsa.numpyTovtkDataArray(tstep.postproc_cell[k], name=f"{k}_cell")
+                )
             for k in postproc_cell_6vector:
-                for i in range(0, 2):
-                    ug.cell_data.add_array(tstep.postproc_cell[k][:, 3*i:3*(i+1)])
-                    ug.cell_data.get_array(counter).name = k + '_' + str(i) + '_cell'
-                    counter += 1
-        ug.cell_data.add_array(coords_a_cell)
-        ug.cell_data.get_array(counter).name = 'coords_a_elem'
-        counter += 1
-
-        ug.point_data.scalars = node_id
-        ug.point_data.scalars.name = 'node_id'
-        point_vector_counter = 1
-        ug.point_data.add_array(local_x, 'vector')
-        ug.point_data.get_array(point_vector_counter).name = 'local_x'
-        point_vector_counter += 1
-        ug.point_data.add_array(local_y, 'vector')
-        ug.point_data.get_array(point_vector_counter).name = 'local_y'
-        point_vector_counter += 1
-        ug.point_data.add_array(local_z, 'vector')
-        ug.point_data.get_array(point_vector_counter).name = 'local_z'
-        point_vector_counter += 1
-        ug.point_data.add_array(coords_a, 'vector')
-        ug.point_data.get_array(point_vector_counter).name = 'coords_a'
+                for i in range(2):
+                    ug.GetCellData().AddArray(
+                        dsa.numpyTovtkDataArray(
+                            tstep.postproc_cell[k][:, 3*i:3*(i+1)], name=f"{k}_{i}_cell"
+                        )
+                    )
+
         if self.settings['include_applied_forces']:
-            point_vector_counter += 1
-            ug.point_data.add_array(app_forces, 'vector')
-            ug.point_data.get_array(point_vector_counter).name = 'app_forces'
-            point_vector_counter += 1
-            ug.point_data.add_array(forces_constraints_nodes, 'vector')
-            ug.point_data.get_array(point_vector_counter).name = 'forces_constraints_nodes'
+            ug.GetPointData().AddArray(dsa.numpyTovtkDataArray(app_forces, name="app_forces"))
+            ug.GetPointData().AddArray(
+                dsa.numpyTovtkDataArray(forces_constraints_nodes, name="forces_constraints_node")
+            )
             if with_gravity:
-                point_vector_counter += 1
-                ug.point_data.add_array(gravity_forces_g[:, 0:3], 'vector')
-                ug.point_data.get_array(point_vector_counter).name = 'gravity_forces'
+                ug.GetPointData().AddArray(
+                    dsa.numpyTovtkDataArray(
+                        gravity_forces_g[:, :3], name="gravity_forces"
+                    )
+                )
 
         if self.settings['include_applied_moments']:
-            point_vector_counter += 1
-            ug.point_data.add_array(app_moment, 'vector')
-            ug.point_data.get_array(point_vector_counter).name = 'app_moments'
-            point_vector_counter += 1
-            ug.point_data.add_array(moments_constraints_nodes, 'vector')
-            ug.point_data.get_array(point_vector_counter).name = 'moments_constraints_nodes'
+            ug.GetPointData().AddArray(
+                dsa.numpyTovtkDataArray(app_moment, name="app_moments")
+            )
+            ug.GetPointData().AddArray(
+                dsa.numpyTovtkDataArray(moments_constraints_nodes, name="moments_constraints_nodes")
+            )
             if with_gravity:
-                point_vector_counter += 1
-                ug.point_data.add_array(gravity_forces_g[:, 3:6], 'vector')
-                ug.point_data.get_array(point_vector_counter).name = 'gravity_moments'
+                ug.GetPointData().AddArray(dsa.numpyTovtkDataArray(gravity_forces_g[:, 3:6], name="gravity_moments"))
+
         if with_postproc_node:
             for k in postproc_node_vector:
-                point_vector_counter += 1
-                ug.point_data.add_array(tstep.postproc_node[k])
-                ug.point_data.get_array(point_vector_counter).name = k + '_point'
+                ug.GetPointData().AddArray(
+                    dsa.numpyTovtkDataArray(
+                        tstep.postproc_node[k], name=f"{k}_point"
+                    )
+                )
             for k in postproc_node_6vector:
-                for i in range(0, 2):
-                    point_vector_counter += 1
-                    ug.point_data.add_array(tstep.postproc_node[k][:, 3*i:3*(i+1)])
-                    ug.point_data.get_array(point_vector_counter).name = k + '_' + str(i) + '_point'
+                for i in range(2):
+                    ug.GetPointData().AddArray(
+                        dsa.numpyTovtkDataArray(
+                            tstep.postproc_node[k][:, 3*i:3*(i+1)], name=f"{k}_{i}_point"
+                        )
+                    )
 
             for k in postproc_node_scalar:
-                point_vector_counter += 1
-                ug.point_data.add_array(tstep.postproc_node[k])
-                ug.point_data.get_array(point_vector_counter).name = k
-
-        write_data(ug, it_filename)
+                ug.GetPointData().AddArray(
+                    dsa.numpyTovtkDataArray(
+                        tstep.postproc_node[k],
+                        name=str(k),
+                    )
+                )
+
+        writer = vtk.vtkXMLPolyDataWriter()
+        writer.SetFileName(it_filename)
+        writer.SetInputData(ug)
+        writer.Write()
 
     def write_for(self, it):
-        it_filename = (self.filename_for +
-                       '%06u' % it)
+        it_filename = f"{self.filename_for}{it:06d}.vtp"
 
-        forces_constraints_FoR = np.zeros((self.data.structure.num_bodies, 3))
-        moments_constraints_FoR = np.zeros((self.data.structure.num_bodies, 3))
-        # TODO: what should I do with the forces of the quaternion?
+        forces_constraints_for = np.zeros((self.data.structure.num_bodies, 3))
+        moments_constraints_for = np.zeros((self.data.structure.num_bodies, 3))
 
         # aero2inertial rotation
         aero2inertial = self.data.structure.timestep_info[it].cga()
 
         # coordinates of corners
-        FoR_coords = np.zeros((self.data.structure.num_bodies, 3))
+        for_coords = np.zeros((self.data.structure.num_bodies, 3))
         if self.settings['include_rbm']:
             offset = np.zeros((3,))
         else:
             offset = self.data.structure.timestep_info[it].mb_FoR_pos[0, 0:3]
-        for ibody in range(self.data.structure.num_bodies):
-            FoR_coords[ibody, :] = self.data.structure.timestep_info[it].mb_FoR_pos[ibody, 0:3] - offset
 
         for ibody in range(self.data.structure.num_bodies):
-            forces_constraints_FoR[ibody, :] = np.dot(aero2inertial,
-                                                  self.data.structure.timestep_info[it].forces_constraints_FoR[ibody, 0:3])
-            moments_constraints_FoR[ibody, :] = np.dot(aero2inertial,
-                                                  self.data.structure.timestep_info[it].forces_constraints_FoR[ibody, 3:6])
-
-        FoRmesh = tvtk.PolyData()
-        FoRmesh.points = FoR_coords
-        for_vector_counter = -1
-        for_vector_counter += 1
-        FoRmesh.point_data.add_array(forces_constraints_FoR , 'vector')
-        FoRmesh.point_data.get_array(for_vector_counter).name = 'forces_constraints_FoR'
-
-        for_vector_counter += 1
-        FoRmesh.point_data.add_array(moments_constraints_FoR , 'vector')
-        FoRmesh.point_data.get_array(for_vector_counter).name = 'moments_constraints_FoR'
-
-        write_data(FoRmesh, it_filename)
+            for_coords[ibody, :] = self.data.structure.timestep_info[it].mb_FoR_pos[ibody, 0:3] - offset
+            forces_constraints_for[ibody, :] = aero2inertial @ self.data.structure.timestep_info[it].forces_constraints_FoR[ibody, :3]
+            moments_constraints_for[ibody, :] = aero2inertial @ self.data.structure.timestep_info[it].forces_constraints_FoR[ibody, 3:6]
+
+        for_mesh = vtk.vtkPolyData(points=for_coords)
+
+        for_mesh.GetPointData().AddArray(
+            dsa.numpyTovtkDataArray(
+                forces_constraints_for,
+                'forces_constraints_FoR',
+            )
+        )
+
+        for_mesh.GetPointData().AddArray(
+            dsa.numpyTovtkDataArray(
+                moments_constraints_for,
+                "moments_constraints_FoR",
+            )
+        )
+
+        writer = vtk.vtkXMLPolyDataWriter()
+        writer.SetFileName(it_filename)
+        writer.SetInputData(for_mesh)
+        writer.Write()
diff --git a/sharpy/rom/utils/krylovutils.py b/sharpy/rom/utils/krylovutils.py
index b6a0b5f89..68490b93a 100644
--- a/sharpy/rom/utils/krylovutils.py
+++ b/sharpy/rom/utils/krylovutils.py
@@ -192,7 +192,7 @@ def construct_krylov(r, lu_A, B, approx_type='Pade', side='b'):
 
 
 def lu_factor(sigma, A):
-    """
+    r"""
     LU Factorisation wrapper of:
 
     .. math:: LU = (\sigma \mathbf{I} - \mathbf{A})
@@ -214,7 +214,7 @@ def lu_factor(sigma, A):
 
 
 def lu_solve(lu_A, b, trans=0):
-    """
+    r"""
     LU solve wrapper.
 
     Computes the solution to
diff --git a/sharpy/sharpy_main.py b/sharpy/sharpy_main.py
index 65f3c0a9e..2d9f0b85c 100644
--- a/sharpy/sharpy_main.py
+++ b/sharpy/sharpy_main.py
@@ -4,7 +4,7 @@
 import warnings
 import sys
 import dill as pickle
-import sharpy.utils.cout_utils as cout
+from sharpy.utils import cout_utils as cout
 from .version import __version__
 
 
diff --git a/sharpy/solvers/modal.py b/sharpy/solvers/modal.py
index 790ce7392..f0cfcc70f 100644
--- a/sharpy/solvers/modal.py
+++ b/sharpy/solvers/modal.py
@@ -209,15 +209,10 @@ def run(self, **kwargs):
 
         # Number of degrees of freedom
         num_str_dof = self.data.structure.num_dof.value
-        if self.rigid_body_motion:
-            num_rigid_dof = 10
-        else:
-            num_rigid_dof = 0
+        num_rigid_dof = 10 if self.rigid_body_motion else 0
 
         num_dof = num_str_dof + num_rigid_dof
 
-        # if NumLambda
-
         # Initialize matrices
         FullMglobal = np.zeros((num_dof, num_dof),
                                dtype=ct.c_double, order='F')
diff --git a/sharpy/solvers/nonlineardynamicmultibody.py b/sharpy/solvers/nonlineardynamicmultibody.py
index bcbd4be7d..d2bec5e6c 100644
--- a/sharpy/solvers/nonlineardynamicmultibody.py
+++ b/sharpy/solvers/nonlineardynamicmultibody.py
@@ -175,7 +175,7 @@ def define_rigid_dofs(self, MB_beam):
             first_dof = last_dof + 0
 
     def assembly_MB_eq_system(self, MB_beam, MB_tstep, ts, dt, Lambda, Lambda_dot, MBdict):
-        """
+        r"""
         This function generates the matrix and vector associated to the linear system to solve a structural iteration
         It usses a Newmark-beta scheme for time integration. Being M, C and K the mass, damping
         and stiffness matrices of the system:
diff --git a/sharpy/structure/utils/modalutils.py b/sharpy/structure/utils/modalutils.py
index 355cc51f2..09dfe7462 100644
--- a/sharpy/structure/utils/modalutils.py
+++ b/sharpy/structure/utils/modalutils.py
@@ -1,8 +1,8 @@
 import numpy as np
 import sharpy.utils.cout_utils as cout
 import sharpy.utils.algebra as algebra
-from tvtk.api import tvtk, write_data
 
+from sharpy.utils.plotutils import plot_frame_to_vtk
 
 def frequency_damping(eigenvalue):
     omega_n = np.abs(eigenvalue)
@@ -41,7 +41,7 @@ def cg(M, use_euler=False):
     return -np.array([Mrr[2, 4], Mrr[0, 5], Mrr[1, 3]]) / Mrr[0, 0]
 
 
-def scale_mode(data, eigenvector, rot_max_deg=15, perc_max=0.15):
+def scale_mode(data, eigenvector, rot_max_deg=15.0, perc_max=0.15):
     """
     Scales the eigenvector such that:
         1) the maximum change in component of the beam cartesian rotation vector
@@ -190,106 +190,48 @@ def write_zeta_vtk(zeta, zeta_ref, filename_root):
     Args:
         zeta (np.array): lattice coordinates to plot
         zeta_ref (np.array): reference lattice used to compute the magnitude of displacements
-        filename_root (str): initial part of filename (full path) without file extension (.vtk)
+        filename_root (str): initial part of filename (full path) without file extension
     """
 
     for i_surf in range(len(zeta)):
 
-        filename = filename_root + "_%02u.vtu" % (i_surf,)
-        _, M, N = zeta[i_surf].shape
-
-        M -= 1
-        N -= 1
-        point_data_dim = (M + 1) * (N + 1)
-        panel_data_dim = M * N
-
-        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)
-        point_struct_id = np.zeros((point_data_dim,), dtype=int)
-        point_struct_mag = np.zeros((point_data_dim,), dtype=float)
-
-        counter = -1
-        # coordinates of corners
-        for i_n in range(N + 1):
-            for i_m in range(M + 1):
-                counter += 1
-                coords[counter, :] = zeta[i_surf][:, i_m, i_n]
-
-        counter = -1
-        node_counter = -1
-        for i_n in range(N + 1):
-            # global_counter = aero.aero2struct_mapping[i_surf][i_n]
-            for i_m in range(M + 1):
-                node_counter += 1
-                # point data
-                # point_struct_id[node_counter]=global_counter
-                point_struct_mag[node_counter] = \
-                    np.linalg.norm(zeta[i_surf][:, i_m, i_n] \
-                                   - zeta_ref[i_surf][:, i_m, i_n])
-
-                if i_n < N and i_m < M:
-                    counter += 1
-                else:
-                    continue
-
-                conn.append([node_counter + 0,
-                             node_counter + 1,
-                             node_counter + M + 2,
-                             node_counter + M + 1])
-                # cell data
-                panel_id[counter] = counter
-                panel_surf_id[counter] = i_surf
-
-        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.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_struct_mag)
-        ug.point_data.get_array(1).name = 'point_displacement_magnitude'
-
-        write_data(ug, filename)
-
-
-def write_modes_vtk(data, eigenvectors, NumLambda, filename_root,
+        filename = f"{filename_root}_{i_surf:02d}.vtu"
+
+        this_zeta = np.swapaxes(zeta[i_surf], 0, -1)
+        this_zeta_ref = np.swapaxes(zeta_ref[i_surf], 0, -1)
+        disp = np.linalg.norm(this_zeta - this_zeta_ref, axis=-1)
+        panel_surf_id = np.ones((this_zeta.shape[0]-1, this_zeta.shape[1]-1), dtype=int) * i_surf
+
+        plot_frame_to_vtk(this_zeta,
+                          filename,
+                          node_scalar_data={'point_displacement_magnitude': disp},
+                            cell_scalar_data={'panel_surface_id': panel_surf_id})
+
+
+def write_modes_vtk(data, eigenvectors, num_lambda, filename_root,
                     rot_max_deg=15., perc_max=0.15, ts=-1):
     """
-    Writes a vtk file for each of the first ``NumLambda`` eigenvectors. When these
+    Writes a vtk file for each of the first ``num_lambda`` eigenvectors. When these
     are associated to the state-space form of the structural equations, only
     the displacement field is saved.
     """
 
-    ### initialise
-    aero = data.aero
+    # initialise
     struct = data.structure
     tsaero = data.aero.timestep_info[ts]
-    tsstr = data.structure.timestep_info[ts]
-
     num_dof = struct.num_dof.value
     eigenvectors = eigenvectors[:num_dof, :]
 
-    # Check whether rigid body motion is selected
-    # Skip rigid body modes
-    if data.settings['Modal']['rigid_body_modes']:
-        num_rigid_body = 10
-    else:
-        num_rigid_body = 0
+    # skip rigid body modes
+    num_rigid_body = 10 if data.settings['Modal']['rigid_body_modes'] else 0
 
-    for mode in range(num_rigid_body, NumLambda):
+    for mode in range(num_rigid_body, num_lambda):
         # scale eigenvector
         eigvec = eigenvectors[:num_dof, mode]
         fact = scale_mode(data, eigvec, rot_max_deg, perc_max)
         eigvec = eigvec * fact
         zeta_mode = get_mode_zeta(data, eigvec)
-        write_zeta_vtk(zeta_mode, tsaero.zeta, filename_root + "_%06u" % (mode,))
+        write_zeta_vtk(zeta_mode, tsaero.zeta, f"{filename_root}_{mode:06d}")
 
 
 def free_modes_principal_axes(phi, mass_matrix, use_euler=False, **kwargs):
diff --git a/sharpy/utils/plotutils.py b/sharpy/utils/plotutils.py
index 7f552f9ab..34d50c532 100644
--- a/sharpy/utils/plotutils.py
+++ b/sharpy/utils/plotutils.py
@@ -1,6 +1,104 @@
 """Plotting utilities
 """
 import numpy as np
+from numpy import ndarray
+
+from typing import Optional
+from os import PathLike
+
+
+import vtk
+from vtk.numpy_interface import algorithms as algs
+from vtk.numpy_interface import dataset_adapter as dsa
+
+def plot_frame_to_vtk(
+    grid_arr: ndarray,
+    filename: str | PathLike,
+    node_scalar_data: Optional[dict[str, ndarray]] = None,
+    node_vector_data: Optional[dict[str, ndarray]] = None,
+    cell_scalar_data: Optional[dict[str, ndarray]] = None,
+    cell_vector_data: Optional[dict[str, ndarray]] = None,
+) -> None:
+    r"""
+    Plot a single timestep of grid data
+    :param grid_arr: Structured grid array with shape (n_x, n_y, n_z, 3)
+    :param filename: Base filename, including directory. Information on the frame number will be
+    appended to this.
+    :param node_scalar_data: Dictionary of node scalar data
+    :param node_vector_data: Dictionary of node vector data
+    :param cell_scalar_data: Dictionary of cell scalar data
+    :param cell_vector_data: Dictionary of cell vector data
+    """
+
+    if grid_arr.shape[-1] != 3:
+        raise ValueError("grid_arr must have trailing dimension of size 3")
+
+    # planar grid should have 3 dimensions, while volume grid should have 4 dimensions
+    match grid_arr.ndim:
+        case 3:
+            is_planar = True
+        case 4:
+            is_planar = False
+        case _:
+            raise ValueError(
+                f"grid_arr must have 3 or 4 dimensions, got {grid_arr.ndim}-D array"
+            )
+
+    sg = vtk.vtkStructuredGrid()
+    if is_planar:
+        sg.SetDimensions(*grid_arr.shape[:-1], 1)
+    else:
+        sg.SetDimensions(*grid_arr.shape[:-1])
+
+    i_swap = 2 - int(is_planar)  # we swap axes as VTK likes z, y, x order
+
+    # add point coordinate data
+    points = vtk.vtkPoints()
+    points_vec = algs.make_vector(
+        *[np.swapaxes(grid_arr[..., i], 0, i_swap).ravel() for i in range(3)]
+    )
+    points.SetData(dsa.numpyTovtkDataArray(points_vec, "Points"))
+    sg.SetPoints(points)
+
+    # cell scalar data
+    if cell_scalar_data is not None:
+        for name, arr in cell_scalar_data.items():
+            sg.GetCellData().AddArray(
+                dsa.numpyTovtkDataArray(np.swapaxes(arr, 0, i_swap).ravel(), name)
+            )
+
+    # cell vector data
+    if cell_vector_data is not None:
+        for name, arr in cell_vector_data.items():
+            vectors = algs.make_vector(
+                np.swapaxes(arr[..., 0], 0, i_swap).ravel(),
+                np.swapaxes(arr[..., 1], 0, i_swap).ravel(),
+                np.swapaxes(arr[..., 2], 0, i_swap).ravel(),
+            )
+            sg.GetCellData().AddArray(dsa.numpyTovtkDataArray(vectors, name))
+
+    # point scalar data
+    if node_scalar_data is not None:
+        for name, arr in node_scalar_data.items():
+            sg.GetPointData().AddArray(
+                dsa.numpyTovtkDataArray(np.swapaxes(arr, 0, i_swap).ravel(), name)
+            )
+
+    # point vector data
+    if node_vector_data is not None:
+        for name, arr in node_vector_data.items():
+            vectors = algs.make_vector(
+                np.swapaxes(arr[..., 0], 0, i_swap).ravel(),
+                np.swapaxes(arr[..., 1], 0, i_swap).ravel(),
+                np.swapaxes(arr[..., 2], 0, i_swap).ravel(),
+            )
+            sg.GetPointData().AddArray(dsa.numpyTovtkDataArray(vectors, name))
+
+    # write to file
+    writer = vtk.vtkXMLStructuredGridWriter()
+    writer.SetFileName(filename)
+    writer.SetInputData(sg)
+    writer.Write()
 
 
 def set_axes_equal(ax):
diff --git a/tests/linear/horten/test_horten.py b/tests/linear/horten/test_horten.py
index 5c8eef3e5..04d106c47 100644
--- a/tests/linear/horten/test_horten.py
+++ b/tests/linear/horten/test_horten.py
@@ -23,10 +23,6 @@ def run_rom_convergence(case_name, case_route='./cases/', output_folder='./outpu
     u_inf = kwargs.get('u_inf', 30)
 
     use_euler = True
-    if use_euler:
-        orient = 'euler'
-    else:
-        orient = 'quat'
 
     case_name += 'M%gN%gMsf%g_u%g' % (M, N, Msf, u_inf)