diff --git a/city_model_structure/attributes/polyhedron.py b/city_model_structure/attributes/polyhedron.py
index 27970eeb..347b5c95 100644
--- a/city_model_structure/attributes/polyhedron.py
+++ b/city_model_structure/attributes/polyhedron.py
@@ -71,6 +71,7 @@ class Polyhedron:
     points_list = surface.points_list
     point_index = 0
     area = surface.area
+    print('AREA:', area)
     while len(triangles) < triangles_count:
       # select a triangle starting at point index
       triangle_points = ' '.join(str(e) for e in [*points_list[point_index:point_index + 9]])
@@ -78,6 +79,8 @@ class Polyhedron:
       rest_points = ' '.join(str(e) for e in [*points_list[0:point_index+3], *points_list[point_index+6:]])
       triangular_surface = Surface(triangle_points, remove_last=False)
       rest_surface = Surface(rest_points, remove_last=False)
+      print('triangular_surface:', triangular_surface.points)
+      print('triangular_surface_AREA:', triangular_surface.area)
       if self._geometry.almost_same_area(area, (triangular_surface.area + rest_surface.area)):
         area = rest_surface.area
         triangles.append(triangular_surface)
diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py
index 44d9bfac..d672691f 100644
--- a/city_model_structure/attributes/surface.py
+++ b/city_model_structure/attributes/surface.py
@@ -222,65 +222,67 @@ class Surface:
       print('NEW METHOD TO CALCULATE AREA')
       print('original:')
       print(self.points)
+      if len(self.points) < 3:
+        area = 0
+      else:
+        # 1. 3D -> 2D
+        z_vector = [0, 0, 1]
+        normal_vector = self.normal
+        print(normal_vector)
+        turning_base_matrix = None
+        points_2d = []
+        x = normal_vector[0]
+        y = normal_vector[1]
+        z = normal_vector[2]
+        print('x:', x)
+        print('y:', y)
+        print('z:', z)
+        if x != 0 or y != 0:
+          # cos(alpha) = n.z/|n|.|z|
+          print('dot_mult:', np.dot(normal_vector, z_vector))
+          alpha = math.acos(np.dot(normal_vector, z_vector) / np.linalg.norm(normal_vector) / np.linalg.norm(z_vector))
+          print('alpha:', alpha)
+          turning_line = np.cross(normal_vector, z_vector)
+          print('turning_line:', turning_line)
+          third_axis = np.cross(normal_vector, turning_line)
+          print('third_axis:', third_axis)
+          # orthonormal base
+          w_1 = turning_line / np.linalg.norm(turning_line)
+          w_2 = normal_vector
+          w_3 = third_axis / np.linalg.norm(third_axis)
+          # turning_base_matrix
+          turning_matrix = np.array([[1, 0, 0],
+                            [0, math.cos(alpha), -math.sin(alpha)],
+                            [0, math.sin(alpha), math.cos(alpha)]])
+          print('turning_matrix:', turning_matrix)
+          base_matrix = np.array([w_1, w_2, w_3])
+          print('base_matrix:', base_matrix)
+          turning_base_matrix = np.matmul(base_matrix.transpose(), turning_matrix.transpose())
+          turning_base_matrix = np.matmul(turning_base_matrix, base_matrix)
+          print('turning_base_matrix:', turning_base_matrix)
 
-      # 1. 3D -> 2D
-      z_vector = [0, 0, 1]
-      normal_vector = self.normal
-      print(normal_vector)
-      turning_base_matrix = None
-      points_2d = []
-      x = normal_vector[0]
-      y = normal_vector[1]
-      z = normal_vector[2]
-      print('x:', x)
-      print('y:', y)
-      print('z:', z)
-      if x != 0 or y != 0:
-        # cos(alpha) = n.z/|n|.|z|
-        print('dot_mult:', np.dot(normal_vector, z_vector))
-        alpha = math.acos(np.dot(normal_vector, z_vector) / np.linalg.norm(normal_vector) / np.linalg.norm(z_vector))
-        print('alpha:', alpha)
-        turning_line = np.cross(normal_vector, z_vector)
-        print('turning_line:', turning_line)
-        third_axis = np.cross(normal_vector, turning_line)
-        print('third_axis:', third_axis)
-        # orthonormal base
-        w_1 = turning_line / np.linalg.norm(turning_line)
-        w_2 = normal_vector
-        w_3 = third_axis / np.linalg.norm(third_axis)
-        # turning_base_matrix
-        turning_matrix = np.array([[1, 0, 0],
-                          [0, math.cos(alpha), -math.sin(alpha)],
-                          [0, math.sin(alpha), math.cos(alpha)]])
-        print('turning_matrix:', turning_matrix)
-        base_matrix = np.array([w_1, w_2, w_3])
-        print('base_matrix:', base_matrix)
-        turning_base_matrix = np.matmul(base_matrix.transpose(), turning_matrix.transpose())
-        turning_base_matrix = np.matmul(turning_base_matrix, base_matrix)
-        print('turning_base_matrix:', turning_base_matrix)
+          if turning_base_matrix is None:
+            print('ERROR')
+          else:
+            for point in self.points:
+              new_point = np.matmul(turning_base_matrix, point)
+              print('new_point:', new_point)
+              points_2d.append(new_point)
+  #            points_2d.append([new_point[0], new_point[1]])
 
-        if turning_base_matrix is None:
-          print('ERROR')
         else:
           for point in self.points:
-            new_point = np.matmul(turning_base_matrix, point)
-            print('new_point:', new_point)
-            points_2d.append(new_point)
-#            points_2d.append([new_point[0], new_point[1]])
+            points_2d.append([point[0], point[1], 0])
 
-      else:
-        for point in self.points:
-          points_2d.append([point[0], point[1], 0])
-
-      polygon_2d = pn.Polygon(np.array(points_2d))
-      print('2D:')
-      print(polygon_2d.points)
-      # 2. calculate area:
-      area = 0
-      for i in range(0, len(polygon_2d.points)-1):
-        point = polygon_2d.points[i]
-        next_point = polygon_2d.points[i+1]
-        area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
+        polygon_2d = pn.Polygon(np.array(points_2d))
+        print('2D:')
+        print(polygon_2d.points)
+        # 2. calculate area:
+        area = 0
+        for i in range(0, len(polygon_2d.points)-1):
+          point = polygon_2d.points[i]
+          next_point = polygon_2d.points[i+1]
+          area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
       self._area = area
     return self._area