diff --git a/CHANGELOG.md b/CHANGELOG.md
index 843ba2cb81a..8f581a15d1b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,6 +19,9 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Changed `compas.scene.Scene` to inherent from `compas.datastructrues.Tree`.
 * Changed `compas.scene.SceneObject` to inherent from `compas.datastructrues.TreeNode`.
 * Changed `compas.geoemetry._core.predicates_3` bug fix in `is_coplanar` while loop when there are 4 points.
+* Changed to implementation of `Mesh.unify_cycles` to use the corresponding function of `compas.topology.orientation`.
+* Fixed bug in `compas.topology.orientation.unify_cycles`.
+* Fixed bug in `Mesh.thickened`.
 
 ### Removed
 
diff --git a/src/compas/datastructures/mesh/mesh.py b/src/compas/datastructures/mesh/mesh.py
index 250e6a8dace..7be50338ef7 100644
--- a/src/compas/datastructures/mesh/mesh.py
+++ b/src/compas/datastructures/mesh/mesh.py
@@ -56,8 +56,8 @@
 from compas.utilities import window
 
 from compas.topology import breadth_first_traverse
-from compas.topology import face_adjacency
 from compas.topology import connected_components
+from compas.topology import unify_cycles
 
 from compas.datastructures.datastructure import Datastructure
 from compas.datastructures.attributes import VertexAttributeView
@@ -2953,42 +2953,25 @@ def unify_cycles(self, root=None):
             The mesh is modified in place.
 
         """
-
-        def unify(node, nbr):
-            # find the common edge
-            for u, v in self.face_halfedges(nbr):
-                if u in self.face[node] and v in self.face[node]:
-                    # node and nbr have edge u-v in common
-                    i = self.face[node].index(u)
-                    j = self.face[node].index(v)
-                    if i == j - 1 or (j == 0 and u == self.face[node][-1]):
-                        # if the traversal of a neighboring halfedge
-                        # is in the same direction
-                        # flip the neighbor
-                        self.face[nbr][:] = self.face[nbr][::-1]
-                        return
-
-        if root is None:
-            root = self.face_sample(size=1)[0]
-
+        vertex_index = {}
+        index_vertex = {}
+        for index, vertex in enumerate(self.vertices()):
+            vertex_index[vertex] = index
+            index_vertex[index] = vertex
         index_face = {index: face for index, face in enumerate(self.faces())}
-        points = self.vertices_attributes("xyz")
-        faces = [self.face_vertices(face) for face in self.faces()]
-
-        adj = face_adjacency(points, faces)
-        adjacency = {}
-        for face in adj:
-            adjacency[index_face[face]] = [index_face[nbr] for nbr in adj[face]]
 
-        visited = breadth_first_traverse(adjacency, root, unify)
+        vertices = self.vertices_attributes("xyz")
+        faces = [[vertex_index[vertex] for vertex in self.face_vertices(face)] for face in self.faces()]
 
-        if len(list(visited)) != self.number_of_faces():
-            raise Exception("Not all faces were visited.")
+        unify_cycles(vertices, faces)
 
         self.halfedge = {key: {} for key in self.vertices()}
-        for fkey in self.faces():
-            for u, v in self.face_halfedges(fkey):
-                self.halfedge[u][v] = fkey
+        for index, vertices in enumerate(faces):
+            face = index_face[index]
+            vertices = [index_vertex[vertex] for vertex in vertices]
+            self.face[face] = vertices
+            for u, v in pairwise(vertices + vertices[:1]):
+                self.halfedge[u][v] = face
                 if u not in self.halfedge[v]:
                     self.halfedge[v][u] = None
 
@@ -4413,6 +4396,11 @@ def face_flatness(self, fkey, maxdev=0.02):
         float
             The flatness.
 
+        Raises
+        ------
+        Exception
+            If the face has more than 4 vertices.
+
         Notes
         -----
         Flatness is computed as the ratio of the distance between the diagonals
@@ -4426,6 +4414,12 @@ def face_flatness(self, fkey, maxdev=0.02):
         """
         vertices = self.face_vertices(fkey)
         f = len(vertices)
+
+        if f == 3:
+            return 0.0
+        if f > 4:
+            raise Exception("Computing face flatness for faces with more than 4 vertices is not supported.")
+
         points = self.vertices_attributes("xyz", keys=vertices) or []
         lengths = [distance_point_point(a, b) for a, b in pairwise(points + points[:1])]
         length = sum(lengths) / f
@@ -5075,8 +5069,8 @@ def thickened(self, thickness=1.0, both=True):
             raise ValueError("Thickness should be a positive number.")
 
         if both:
-            mesh_top = self.offset(+0.5 * thickness)
-            mesh_bottom = self.offset(-0.5 * thickness)
+            mesh_top = self.offset(+0.5 * thickness)  # type: Mesh
+            mesh_bottom = self.offset(-0.5 * thickness)  # type: Mesh
         else:
             mesh_top = self.offset(thickness)
             mesh_bottom = self.copy()
@@ -5085,7 +5079,8 @@ def thickened(self, thickness=1.0, both=True):
         mesh_bottom.flip_cycles()
 
         # join parts
-        thickened_mesh = mesh_top.join(mesh_bottom)
+        thickened_mesh = mesh_top.copy()  # type: Mesh
+        thickened_mesh.join(mesh_bottom)
 
         # close boundaries
         n = thickened_mesh.number_of_vertices() / 2
diff --git a/src/compas/topology/orientation.py b/src/compas/topology/orientation.py
index 4593a58e816..20a1f3c8712 100644
--- a/src/compas/topology/orientation.py
+++ b/src/compas/topology/orientation.py
@@ -88,7 +88,7 @@ def _face_adjacency(vertices, faces, nmax=10, radius=10.0):
                         found.add(nbr)
                         break
 
-        adjacency[face] = nbrs
+        adjacency[index] = nbrs
 
     return adjacency
 
@@ -98,8 +98,10 @@ def face_adjacency(points, faces):
 
     Parameters
     ----------
-    mesh : :class:`compas.datastructures.Mesh`
-        A mesh object.
+    points : list[point]
+        The vertex locations of the faces.
+    faces : list[list[int]]
+        The faces defined as list of indices in the points list.
 
     Returns
     -------