diff --git a/hypervehicle/components/component.py b/hypervehicle/components/component.py
index f432559..3d0522f 100644
--- a/hypervehicle/components/component.py
+++ b/hypervehicle/components/component.py
@@ -1,21 +1,20 @@
-import os
-import time
 import pymeshfix
 import numpy as np
 from stl import mesh
 import multiprocess as mp
 from copy import deepcopy
 from abc import abstractmethod
+from collections import Counter
 from hypervehicle.geometry import Vector3
 from gdtk.geom.sgrid import StructuredGrid
 from typing import Callable, Union, Optional
+from hypervehicle.utilities import surfce_to_stl
 from hypervehicle.geometry import (
     CurvedPatch,
     RotatedPatch,
     MirroredPatch,
     OffsetPatchFunction,
 )
-from hypervehicle.utilities import parametricSurfce2stl
 
 
 class AbstractComponent:
@@ -90,6 +89,7 @@ class Component(AbstractComponent):
     def __init__(
         self,
         params: dict = None,
+        edges: list = None,
         stl_resolution: int = 2,
         verbosity: int = 1,
         name: str = None,
@@ -100,8 +100,13 @@ def __init__(
         # Save parameters
         self.params = params
 
+        # Save edges
+        self.edges = edges if edges is not None else []
+
         # Processed objects
         self.patches = {}  # Parametric patches (continuous)
+        self.patch_res_r = {}  # corresponding stl resolution in 'r' direction
+        self.patch_res_s = {}  # corresponding stl resolution in 's' direction
 
         # VTK Attributes
         self.grids = None  # Structured grids
@@ -252,32 +257,44 @@ def surface(self, resolution: int = None):
                 "No patches have been generated. " + "Please call .generate_patches()."
             )
 
-        # Prepare multiprocessing arguments iterable
-        def wrapper(key: str, patch):
-            flip = True if key.split("_")[-1] == "mirrored" else False
-            res = stl_resolution
-
-            if "swept" in key:
-                # Swept fuselage component
-                res = (
-                    int(stl_resolution / 4)
-                    if "end" in key
-                    else int(stl_resolution / 4) * 4
-                )
-                flip = True if "1" in key else False
+        # Create case list
+        if isinstance(stl_resolution, int):
+            case_list = [
+                [k, patch, stl_resolution, stl_resolution]
+                for k, patch in self.patches.items()
+            ]
+        else:
+            case_list = [
+                [k, patch, self.patch_res_r[k], self.patch_res_r[k]]
+                for k, patch in self.patches.items()
+            ]
 
-            surface = parametricSurfce2stl(
-                patch, res, flip_faces=flip, **self._clustering
-            )
-            return (key, surface)
+        # Prepare multiprocessing arguments iterable
+        def wrapper(key: str, patch, res_r: int, res_s: int):
+            surface = surfce_to_stl(patch, res_r, res_s, **self._clustering)
+            return surface
 
-        # Initialise surfaces and pool
         self.surfaces = {}
-        pool = mp.Pool()
+        multiprocess = False  # flag to disable multiprocessing for debugging
+        # TODO - move multiprocess to arg / config option
+        if multiprocess is True:
+            # Initialise surfaces and pool
+            pool = mp.Pool()
+
+            # Submit tasks
+            print(f"START: Creating stl - multiprocessor run.")
+            for result in pool.starmap(wrapper, case_list):
+                self.surfaces[result[0]] = result[1]
+            print("  DONE: Creating stl - multiprocess.")
 
-        # Submit tasks
-        for result in pool.starmap(wrapper, self.patches.items()):
-            self.surfaces[result[0]] = result[1]
+        else:
+            for case in case_list:
+                k = case[0]
+                pat = case[1]
+                print(f"START: Creating stl for '{k}'.")
+                result = wrapper(k, pat, case[2], case[3])
+                self.surfaces[result[0]] = result[1]
+                print("  DONE: Creating stl.")
 
     def to_vtk(self):
         raise NotImplementedError("This method has not been implemented yet.")
@@ -332,3 +349,140 @@ def add_clustering_options(
 
         if j_clustering_func:
             self._clustering.update({"j_clustering_func": j_clustering_func})
+
+    def stl_check(
+        self, project_area: Optional[bool] = True, matching_lines: Optional[bool] = True
+    ):
+        """Check the STL mesh."""
+        # TODO - add comment annotations
+        pass_flag = True
+        mesh = self.mesh
+        print(f"    MESH CHECK")
+        print(f"        mesh: {mesh}")
+        if project_area:
+            print(f"    Project Area Check")
+            normals = np.asarray(mesh.normals, dtype=np.float64)
+            allowed_max_errors = np.abs(normals).sum(axis=0) * np.finfo(np.float32).eps
+            sum_normals = normals.sum(axis=0)
+            print(f"        allowed error: {allowed_max_errors}")
+            print(f"        normals sum: {sum_normals}")
+            if ((np.abs(sum_normals)) <= allowed_max_errors).all():
+                print(f"      PASS")
+            else:
+                print(f"      FAIL")
+                pass_flag = False
+
+        if matching_lines:
+            print(f"    Matching Lines Check")
+            reversed_triangles = (
+                np.cross(mesh.v1 - mesh.v0, mesh.v2 - mesh.v0) * mesh.normals
+            ).sum(axis=1) < 0
+            import itertools
+
+            directed_edges = {
+                tuple(edge.ravel() if not rev else edge[::-1, :].ravel())
+                for rev, edge in zip(
+                    itertools.cycle(reversed_triangles),
+                    itertools.chain(
+                        mesh.vectors[:, (0, 1), :],
+                        mesh.vectors[:, (1, 2), :],
+                        mesh.vectors[:, (2, 0), :],
+                    ),
+                )
+            }
+        undirected_edges = {frozenset((edge[:3], edge[3:])) for edge in directed_edges}
+        edge_check = len(directed_edges) == 3 * mesh.data.size
+        print(f"        len(directed_edges) == 3 * mesh.data.size")
+        if edge_check:
+            print(f"        {len(directed_edges)} == {3*mesh.data.size}")
+            print(f"      PASS")
+        else:
+            print(f"        {len(directed_edges)} != {3*mesh.data.size}")
+            print(f"      FAIL")
+            print(f"        The number of edges should be N_cells*3")
+            print(f"        len(directed_edges)={len(directed_edges)}")
+            print(f"        mesh.data.size={mesh.data.size}")
+            pass_flag = False
+        pair_check = len(directed_edges) == 2 * len(undirected_edges)
+        print(f"        len(directed_edges) == 2 * len(undirected_edges)")
+        if pair_check:
+            print(f"        {len(directed_edges)} == {2 * len(undirected_edges)}")
+            print(f"      PASS")
+            return pass_flag
+        else:
+            print(f"        {len(directed_edges)} != {2 * len(undirected_edges)}")
+            print(f"      FAIL")
+            print(f"        All edges should be in pairs")
+            print(f"        len_directed:{len(directed_edges)}")
+            print(f"        len_undirect (pairs+leftover):{len(undirected_edges)}")
+            pass_flag = False
+
+        edge_list = []
+        for edge in directed_edges:
+            edge_list.append(frozenset((edge[:3], edge[3:])))
+
+        counter_out = Counter(edge_list)
+
+        k_list = []
+        for k, v in counter_out.items():
+            if v == 2:
+                k_list.append(k)
+        for k in k_list:
+            del counter_out[k]
+
+        return pass_flag
+
+    def stl_repair(self, small_distance: float = 1e-6):
+        """Attempts to repair stl mesh issues.
+
+        Parameters
+        -----------
+        small_distance : Float, optional
+            Vectors less than this distance apart will be set to the same value.
+        """
+        mesh = self.mesh
+        N_faces = np.shape(mesh)[0]
+        vectors = np.vstack((mesh.v0, mesh.v1, mesh.v2))
+
+        def find_groups(vector, small_distance):
+            """find indices for groups of points that within small_distance from one another."""
+            sort_i = np.argsort(vector)
+            groups = []
+            li = []
+            count = 0
+            for i in range(len(vector) - 1):
+                if count == 0:
+                    x0 = vector[sort_i[i]]
+                x1 = vector[sort_i[i + 1]]
+                if abs(x1 - x0) < small_distance:
+                    if count == 0:
+                        li.append(sort_i[i])
+                        li.append(sort_i[i + 1])
+                        count = 1
+                    else:
+                        li.append(sort_i[i + 1])
+                else:
+                    groups.append(li)
+                    li = []
+                    count = 0
+                if i == len(vector) - 2 and count > 0:
+                    groups.append(li)
+            return groups
+
+        iix = find_groups(vectors[:, 0], small_distance)
+        iiy = find_groups(vectors[:, 1], small_distance)
+        iiz = find_groups(vectors[:, 2], small_distance)
+
+        # find intersecting sets and take average
+        for ix in iix:
+            for iy in iiy:
+                common0 = set(ix) & set(iy)
+                if len(common0) == 0:
+                    continue  # jumpt to next loop iteration if only one entry
+                for iz in iiz:
+                    common = list(common0 & set(iz))
+                    if common:
+                        vectors[common] = np.mean(vectors[common], axis=0)
+        mesh.v0 = vectors[0:N_faces, :]
+        mesh.v1 = vectors[N_faces : 2 * N_faces, :]
+        mesh.v2 = vectors[2 * N_faces : 3 * N_faces, :]
diff --git a/hypervehicle/components/fin.py b/hypervehicle/components/fin.py
index 3c101ae..a94edbd 100644
--- a/hypervehicle/components/fin.py
+++ b/hypervehicle/components/fin.py
@@ -4,7 +4,7 @@
 from gdtk.geom.vector3 import Vector3
 from typing import Callable, Optional
 from gdtk.geom.surface import CoonsPatch
-from gdtk.geom.path import Line, Polyline
+from gdtk.geom.path import Line, Polyline, ReversedPath
 from hypervehicle.components.component import Component
 from hypervehicle.components.constants import FIN_COMPONENT
 from hypervehicle.geometry import (
@@ -48,6 +48,16 @@ def __init__(
     ) -> None:
         """Creates a new fin component.
 
+        Fin geometry defined by 4 points and straight edges
+        between the points that define the fin planform.
+        Leading Edge runs p3->p2->p1.
+
+        p1--N--p2
+        |         \
+        w           e      <---- FLOW
+        |             \
+        p0-----S------p3
+
         Parameters
         ----------
         p0 : Vector3
@@ -135,7 +145,9 @@ def __init__(
             "offset_function": offset_func,
         }
 
-        super().__init__(params, stl_resolution, verbosity, name)
+        super().__init__(
+            params=params, stl_resolution=stl_resolution, verbosity=verbosity, name=name
+        )
 
     def generate_patches(self):
         # Initialise
@@ -280,11 +292,17 @@ def generate_patches(self):
 
         p0p0_top = Line(p0=p0, p1=p0 - Vector3(0, 0, fin_thickness / 2))
         p0_botp0 = Line(p0=p0 + Vector3(0, 0, fin_thickness / 2), p1=p0)
-        p3p3_top = Line(p0=p3, p1=p3 - Vector3(0, 0, fin_thickness / 2))
-        p3_botp3 = Line(p0=p3 + Vector3(0, 0, fin_thickness / 2), p1=p3)
+        # p3p3_top_reversed = SubRangedPath(p3p3_top, 1, 0)
+        # p3p3_bot_reversed = SubRangedPath(p3p3_bot, 1, 0)
+        p3p3_top_reversed = ReversedPath(p3p3_top)
+        p3p3_bot_reversed = ReversedPath(p3p3_bot)
 
-        bot_1 = CoonsPatch(north=p3p0_top, south=p3p0, east=p0p0_top, west=p3p3_top)
-        bot_2 = CoonsPatch(north=p3p0, south=p3p0_bot, east=p0_botp0, west=p3_botp3)
+        bot_1 = CoonsPatch(
+            north=p3p0_top, south=p3p0, east=p0p0_top, west=p3p3_top_reversed
+        )
+        bot_2 = CoonsPatch(
+            north=p3p0, south=p3p0_bot, east=p0_botp0, west=p3p3_bot_reversed
+        )
 
         temp_fin_patch_dict["bot_ellip"] = bottom_ellipse_patch
         temp_fin_patch_dict["bot_1"] = bot_1
diff --git a/hypervehicle/components/swept.py b/hypervehicle/components/swept.py
index 4416269..51bd2c5 100644
--- a/hypervehicle/components/swept.py
+++ b/hypervehicle/components/swept.py
@@ -1,6 +1,8 @@
-from typing import List
+from typing import List, Dict, Optional
+from gdtk.geom.path import ReversedPath
+from gdtk.geom.surface import CoonsPatch
+from hypervehicle.geometry import SweptPatchfromEdges
 from hypervehicle.components.component import Component
-from hypervehicle.geometry import SweptPatch, CoonsPatch
 from hypervehicle.components.constants import SWEPT_COMPONENT
 
 
@@ -9,32 +11,122 @@ class SweptComponent(Component):
 
     def __init__(
         self,
-        cross_sections: List[CoonsPatch],
-        sweep_axis: str = "z",
-        stl_resolution: int = 2,
-        verbosity: int = 1,
-        name: str = None,
+        cross_sections: List[List],
+        close_ends: Optional[bool] = True,
+        stl_resolution: int | Dict[str, int] = 1,
+        verbosity: Optional[int] = 1,
+        name: Optional[str] = None,
     ) -> None:
         """Create a swept component.
 
         Parameters
         ----------
-        cross_sections : list, optional
-            A list of cross-sectional patches to sweep through.
-        sweep_axis : str, optional
-            The axis to sweep the cross sections through. The default
-            is z.
+        cross_sections : list
+            A list containing cross-sections. Each cross-section is a list of
+            paths that define one cross-section.
+
+        close_ends : bool, optional
+            If true the first and last cross-section will be used to close
+            the sweap component. Only supported is cross-section is defined by
+            4 paths. The default is True.
+
+        stl_resolution : int | dict[str, int], optional
+            Defines different stl resolution for different edges.  Keys are
+            are: 'e0', 'e1', ... 'eN', 'sweep' for edges in first cross-section
+            and for swept edges. The default is 1.
         """
-        self.cross_sections = cross_sections
-        self.sweep_axis = sweep_axis
         super().__init__(stl_resolution=stl_resolution, verbosity=verbosity, name=name)
+        self.cross_sections = cross_sections
+        self.close_ends = close_ends
+        self.n_slices = len(self.cross_sections)
+        self.n_edges = len(self.cross_sections[0])
+        self._check_options()
+
+    def _check_options(self, small_number: Optional[int] = 1e-12):
+        for ns, cs in enumerate(self.cross_sections):
+            if len(cs) != self.n_edges:
+                # check that each c/s has correct number of edges
+                raise Exception(
+                    f"Swept Component {self.name}."
+                    + f"Slice {ns} has incorrect number of edges.\n"
+                    + f"N_e={len(cs)} - {self.n_edges} expected."
+                )
+            for i in range(len(cs)):
+                # check that edges in each c/s form a closed loop
+                if i < len(cs) - 1:
+                    ip = i + 1
+                else:
+                    ip = 0
+                p1 = cs[i](1)
+                p0 = cs[ip](0)
+                if (
+                    abs(p0.x - p1.x) > small_number
+                    or abs(p0.y - p1.y) > small_number
+                    or abs(p0.z - p1.z) > small_number
+                ):
+                    raise Exception(
+                        f"Swept Component {self.name}, Slice {ns}, edges not closed.\n"
+                        + f"edges[{i}](1) != edges[{ip}](0)\n"
+                        + f"{p1} != {p0}"
+                    )
+        if self.close_ends is True and self.n_edges != 4:
+            raise Exception(
+                f"Swept Component {self.name}. Combination of "
+                + f"close_ends={self.close_ends} and N_edge={self.n_edges} is "
+                + f"not supported."
+            )
+        if self.close_ends and isinstance(self.stl_resolution, Dict):
+            flag = 0
+            if self.stl_resolution["e0"] != self.stl_resolution["e2"]:
+                print("edge 'e0' and 'e2' don't have same stl_resolution.")
+                flag = 1
+            if self.stl_resolution["e1"] != self.stl_resolution["e3"]:
+                print("edge 'e1' and 'e3' don't have same stl_resolution.")
+                flag = 1
+            if flag > 0:
+                raise Exception(f"stl_resolution not compatible for close_end.")
 
     def generate_patches(self):
-        p = SweptPatch(
-            cross_sections=self.cross_sections,
-            sweep_axis=self.sweep_axis,
-        )
-
-        self.patches["swept_patch"] = p
-        self.patches["swept_patch_end_0"] = self.cross_sections[0]
-        self.patches["swept_patch_end_1"] = self.cross_sections[-1]
+        for ne in range(self.n_edges):
+            k = f"swept_patch_{ne}"
+            edges = []
+            for cs in self.cross_sections:
+                edges.append(cs[ne])
+            p = SweptPatchfromEdges(edges=edges)
+            self.patches[k] = p
+            if isinstance(self.stl_resolution, int):
+                self.patch_res_r[k] = self.stl_resolution
+                self.patch_res_s[k] = self.stl_resolution
+            else:
+                self.patch_res_r[k] = self.stl_resolution["sweep"]
+                self.patch_res_s[k] = self.stl_resolution[f"e{ne}"]
+
+        if self.close_ends == True:
+            edges = self.cross_sections[0]  # front
+            south = edges[0]
+            east = edges[1]
+            north = ReversedPath(edges[2])
+            west = ReversedPath(edges[3])
+            self.patches["swept_patch_end_0"] = CoonsPatch(
+                south=south, north=north, west=west, east=east
+            )
+            if isinstance(self.stl_resolution, int):
+                self.patch_res_r["swept_patch_end_0"] = self.stl_resolution
+                self.patch_res_s["swept_patch_end_0"] = self.stl_resolution
+            else:
+                self.patch_res_r["swept_patch_end_0"] = self.stl_resolution["e0"]
+                self.patch_res_s["swept_patch_end_0"] = self.stl_resolution["e1"]
+            edges = self.cross_sections[-1]  # rear
+            south = ReversedPath(edges[0])
+            east = ReversedPath(edges[3])
+            north = edges[2]
+            west = edges[1]
+            self.patches["swept_patch_end_1"] = CoonsPatch(
+                south=south, north=north, west=west, east=east
+            )
+            if isinstance(self.stl_resolution, int):
+                self.patch_res_r["swept_patch_end_1"] = self.stl_resolution
+                self.patch_res_s["swept_patch_end_1"] = self.stl_resolution
+            else:
+                self.patch_res_r["swept_patch_end_1"] = self.stl_resolution["e0"]
+                self.patch_res_s["swept_patch_end_1"] = self.stl_resolution["e1"]
diff --git a/hypervehicle/components/wing.py b/hypervehicle/components/wing.py
index 4d9d709..be01c60 100644
--- a/hypervehicle/components/wing.py
+++ b/hypervehicle/components/wing.py
@@ -151,7 +151,9 @@ def __init__(
             "CLOSE_WING": close_wing,
         }
 
-        super().__init__(params, stl_resolution, verbosity, name)
+        super().__init__(
+            params=params, stl_resolution=stl_resolution, verbosity=verbosity, name=name
+        )
 
         # Extract construction points for planform
         # TODO - avoid pre-defined params dict structure for flexibility
diff --git a/hypervehicle/geometry.py b/hypervehicle/geometry.py
index 1a0114e..2df37d6 100644
--- a/hypervehicle/geometry.py
+++ b/hypervehicle/geometry.py
@@ -682,6 +682,44 @@ def __call__(self, r, s) -> Vector3:
         return Vector3(**args)
 
 
+class SweptPatchfromEdges(ParametricSurface):
+    """Creates a swept patch from a series of cross sections.
+    Cross sections do not need to be parallel or aligned."""
+
+    __slots__ = ["edges"]
+
+    def __init__(self, edges: list) -> None:
+        """Construct the SweptPatch object.
+
+        Parameters
+        -----------
+        edges : list
+            A list containing the edges through which to sweep.
+        """
+        self.edges = edges
+        self.n_edges = len(edges)
+
+    def __repr__(self):
+        return "Swept Patch"
+
+    def __call__(self, r, s) -> Vector3:
+        r_slices = np.linspace(0, 1, self.n_edges)
+        if r in r_slices:
+            # can evaluate exactly
+            i = np.where(r_slices == r)[0][0]
+            return self.edges[i](s)
+        else:
+            # need to interpolate
+            index_above = np.where(r_slices > r)[0][0]
+            index_below = index_above - 1
+            above = self.edges[index_above](s)
+            below = self.edges[index_below](s)
+            alpha = (r - r_slices[index_below]) / (
+                r_slices[index_above] - r_slices[index_below]
+            )
+            return (1 - alpha) * below + alpha * above
+
+
 class RotatedPatch(ParametricSurface):
     """
     Rotates a surface about a point in an axis-specified direction.
diff --git a/hypervehicle/utilities.py b/hypervehicle/utilities.py
index e7d42cf..4590c8d 100644
--- a/hypervehicle/utilities.py
+++ b/hypervehicle/utilities.py
@@ -11,15 +11,16 @@
 from typing import Dict, List, Optional
 
 
-def parametricSurfce2stl(
+def surfce_to_stl(
     parametric_surface,
-    triangles_per_edge: int,
+    triangles_per_edge_r: int,
+    triangles_per_edge_s: int,
     si: float = 1.0,
     sj: float = 1.0,
     mirror_y=False,
-    flip_faces=False,
     i_clustering_func: callable = None,
     j_clustering_func: callable = None,
+    verbosity=0,
 ) -> mesh.Mesh:
     """
     Function to convert parametric_surface generated using the Eilmer Geometry
@@ -36,8 +37,11 @@ def parametricSurfce2stl(
         sj : float, optional
             The clustering in the j-direction. The default is 1.0.
 
-        triangles_per_edge : int
-            The resolution for the stl object.
+        triangles_per_edge_r : int
+            The resolution for the stl object - x/i direction.
+
+        triangles_per_edge_s : int
+            The resolution for the stl object - y/j direction.
 
         mirror_y : bool, optional
             Create mirror image about x-z plane. The default is False.
@@ -55,10 +59,11 @@ def parametricSurfce2stl(
     stl_mesh : Mesh
         The numpy-stl mesh.
     """
-    # TODO - allow different ni and nj discretisation
+    ni = triangles_per_edge_r
+    nj = triangles_per_edge_s
 
-    ni = triangles_per_edge
-    nj = triangles_per_edge
+    if verbosity > 1:
+        print(f"Triangles per edge: (ni, nj) = ({ni}, {nj})")
 
     # Create list of vertices
     if i_clustering_func:
@@ -73,83 +78,57 @@ def parametricSurfce2stl(
 
     y_mult: int = -1 if mirror_y else 1
 
-    # Create vertices for corner points of each quad cell
-    # columns x, y, z for each vertex row
-    # quad exterior vertices + quad centres
-    vertices: np.ndarray = np.zeros(((ni + 1) * (nj + 1) + ni * nj, 3))
-    centre_ix: int = (ni + 1) * (nj + 1)
-
+    if verbosity > 1:
+        print(f"r_list = {r_list}")
+        print(f"s_list = {s_list}")
+
+    # For stl generation we split each cell into 4x triangles with indece [0, 3]
+    #   p01-------p11
+    #    | \  2  / |
+    #    |   \  /  |
+    #    | 3  pc 1 |
+    #    |   /  \  |
+    #    |  /  0  \ |
+    #   p00-------p10
+
+    N_triangles = 4 * (ni) * (nj)
+    stl_mesh = mesh.Mesh(np.zeros(N_triangles, dtype=mesh.Mesh.dtype))
+    # mesh.vectors contains a list defining three corners of each triangle.
+    t = 0
     # For vertices along the x direction (i)
-    for i, r in enumerate(r_list):
-        # For vertices along the y direction (j)
-        for j, s in enumerate(s_list):
-            # Evaluate position
-            pos = parametric_surface(r, s)
-
-            # Assign vertex
-            vertices[j * (ni + 1) + i] = np.array([pos.x, y_mult * pos.y, pos.z])
-
-            # Create vertices for centre point of each quad cell
-            try:
-                # Try index to bail before calling surface
-                vertices[centre_ix + (j * ni + i)]
-
-                r0 = r_list[i]
-                r1 = r_list[i + 1]
-                s0 = s_list[j]
-                s1 = s_list[j + 1]
-
-                # Get corner points
-                pos00 = parametric_surface(r0, s0)
-                pos10 = parametric_surface(r1, s0)
-                pos01 = parametric_surface(r0, s1)
-                pos11 = parametric_surface(r1, s1)
-
-                # Evaluate quad centre coordinate
-                pos_x = 0.25 * (pos00.x + pos10.x + pos01.x + pos11.x)
-                pos_y = 0.25 * (pos00.y + pos10.y + pos01.y + pos11.y)
-                pos_z = 0.25 * (pos00.z + pos10.z + pos01.z + pos11.z)
-
-                # Assign quad centre vertices
-                vc = np.array([pos_x, y_mult * pos_y, pos_z])
-                vertices[centre_ix + (j * ni + i)] = vc
-
-            except IndexError:
-                # Index out of bounds
-                pass
-
-    # Create list of faces, defining the face vertices
-    faces = []
     for i in range(ni):
+        # For vertices along the y direction (j)
         for j in range(nj):
-            p00 = j * (nj + 1) + i  # bottom left
-            p10 = j * (nj + 1) + i + 1  # bottom right
-            p01 = (j + 1) * (ni + 1) + i  # top left
-            p11 = (j + 1) * (ni + 1) + i + 1  # top right
-
-            pc = centre_ix + j * min(ni, nj) + i  # vertex at centre of cell
-
-            if mirror_y or flip_faces:
-                faces.append([p00, pc, p10])
-                faces.append([p10, pc, p11])
-                faces.append([p11, pc, p01])
-                faces.append([p01, pc, p00])
-            else:
-                faces.append([p00, p10, pc])
-                faces.append([p10, p11, pc])
-                faces.append([p11, p01, pc])
-                faces.append([p01, p00, pc])
-
-    faces = np.array(faces)
-
-    # Create the STL mesh object
-    stl_mesh = mesh.Mesh(np.zeros(faces.shape[0], dtype=mesh.Mesh.dtype))
-    for ix, face in enumerate(faces):
-        # For each face
-        for c in range(3):
-            # For each coordinate x/y/z
-            stl_mesh.vectors[ix][c] = vertices[face[c], :]
-
+            # get corner points for current cell (quad)
+            pos00 = parametric_surface(r_list[i], s_list[j])
+            pos10 = parametric_surface(r_list[i + 1], s_list[j])
+            pos01 = parametric_surface(r_list[i], s_list[j + 1])
+            pos11 = parametric_surface(r_list[i + 1], s_list[j + 1])
+            # get centre for current cell (quad)
+            pos_x = 0.25 * (pos00.x + pos10.x + pos01.x + pos11.x)
+            pos_y = 0.25 * (pos00.y + pos10.y + pos01.y + pos11.y)
+            pos_z = 0.25 * (pos00.z + pos10.z + pos01.z + pos11.z)
+
+            # add triangle 0 [p00, p10, pc]
+            stl_mesh.vectors[t][0] = np.array([pos00.x, y_mult * pos00.y, pos00.z])
+            stl_mesh.vectors[t][1] = np.array([pos10.x, y_mult * pos10.y, pos10.z])
+            stl_mesh.vectors[t][2] = np.array([pos_x, y_mult * pos_y, pos_z])
+            t += 1
+            # add triangle 1 [p10, p11, pc]
+            stl_mesh.vectors[t][0] = np.array([pos10.x, y_mult * pos10.y, pos10.z])
+            stl_mesh.vectors[t][1] = np.array([pos11.x, y_mult * pos11.y, pos11.z])
+            stl_mesh.vectors[t][2] = np.array([pos_x, y_mult * pos_y, pos_z])
+            t += 1
+            # add triangle 2 [p11, p01, pc]
+            stl_mesh.vectors[t][0] = np.array([pos11.x, y_mult * pos11.y, pos11.z])
+            stl_mesh.vectors[t][1] = np.array([pos01.x, y_mult * pos01.y, pos01.z])
+            stl_mesh.vectors[t][2] = np.array([pos_x, y_mult * pos_y, pos_z])
+            t += 1
+            # add triangle 3 [p01, p00, pc]
+            stl_mesh.vectors[t][0] = np.array([pos01.x, y_mult * pos01.y, pos01.z])
+            stl_mesh.vectors[t][1] = np.array([pos00.x, y_mult * pos00.y, pos00.z])
+            stl_mesh.vectors[t][2] = np.array([pos_x, y_mult * pos_y, pos_z])
+            t += 1
     return stl_mesh
 
 
diff --git a/tests/test_sweep.py b/tests/test_sweep.py
index 818e687..5379fc3 100644
--- a/tests/test_sweep.py
+++ b/tests/test_sweep.py
@@ -1,7 +1,7 @@
 import os
 import numpy as np
 from stl import Mesh
-from hypervehicle.utilities import parametricSurfce2stl
+from hypervehicle.utilities import surfce_to_stl
 from hypervehicle.geometry import SweptPatch, CoonsPatch, Vector3
 
 
@@ -26,9 +26,9 @@ def test_sweep():
 
     p = SweptPatch(sections)
 
-    s = parametricSurfce2stl(p, 20)
-    s1 = parametricSurfce2stl(c1, 20)
-    s2 = parametricSurfce2stl(c2, 20)
+    s = surfce_to_stl(p, 20)
+    s1 = surfce_to_stl(c1, 20)
+    s2 = surfce_to_stl(c2, 20)
 
     stl_data = [s.data, s1.data, s2.data]
     m = Mesh(np.concatenate(stl_data))