diff --git a/hub/city_model_structure/attributes/polygon.py b/hub/city_model_structure/attributes/polygon.py
index 94d47725..2755d399 100644
--- a/hub/city_model_structure/attributes/polygon.py
+++ b/hub/city_model_structure/attributes/polygon.py
@@ -9,9 +9,13 @@ from __future__ import annotations
 import math
 import sys
 from typing import List
+from hub.hub_logger import logger
 import numpy as np
 from trimesh import Trimesh
 import trimesh.intersections
+import trimesh.creation
+import trimesh.geometry
+from shapely.geometry.polygon import Polygon as shapley_polygon
 
 from hub.city_model_structure.attributes.plane import Plane
 from hub.city_model_structure.attributes.point import Point
@@ -22,6 +26,7 @@ class Polygon:
   """
   Polygon class
   """
+  # todo: review with @Guille: Points, Coordinates, Vertices, Faces
 
   def __init__(self, coordinates):
     self._area = None
@@ -66,20 +71,6 @@ class Polygon:
     """
     return self._coordinates
 
-  @staticmethod
-  def _module(vector):
-    x2 = vector[0] ** 2
-    y2 = vector[1] ** 2
-    z2 = vector[2] ** 2
-    return math.sqrt(x2+y2+z2)
-
-  @staticmethod
-  def _scalar_product(vector_0, vector_1):
-    x = vector_0[0] * vector_1[0]
-    y = vector_0[1] * vector_1[1]
-    z = vector_0[2] * vector_1[2]
-    return x+y+z
-
   def contains_point(self, point):
     """
     Determines if the given point is contained by the current polygon
@@ -98,9 +89,9 @@ class Polygon:
       vector_1[0] = vector_1[0] - point.coordinates[0]
       vector_1[1] = vector_1[1] - point.coordinates[1]
       vector_1[2] = vector_1[2] - point.coordinates[2]
-      module = Polygon._module(vector_0) * Polygon._module(vector_1)
+      module = np.linalg.norm(vector_0) * np.linalg.norm(vector_1)
 
-      scalar_product = Polygon._scalar_product(vector_0, vector_1)
+      scalar_product = np.dot(vector_0, vector_1)
       angle = np.pi/2
       if module != 0:
         angle = abs(np.arcsin(scalar_product / module))
@@ -150,69 +141,17 @@ class Polygon:
     Get surface area in square meters
     :return: float
     """
-    # New method to calculate area
     if self._area is None:
-      if len(self.points) < 3:
-        sys.stderr.write('Warning: the area of a line or point cannot be calculated 1. Area = 0\n')
-        return 0
-      alpha = 0
-      vec_1 = self.points[1].coordinates - self.points[0].coordinates
-      for i in range(2, len(self.points)):
-        vec_2 = self.points[i].coordinates - self.points[0].coordinates
-        alpha += self._angle_between_vectors(vec_1, vec_2)
-      if alpha == 0:
-        sys.stderr.write('Warning: the area of a line or point cannot be calculated 2. Area = 0\n')
-        return 0
-      horizontal_points = self._points_rotated_to_horizontal
-      area = 0
-      for i in range(0, len(horizontal_points) - 1):
-        point = horizontal_points[i]
-        next_point = horizontal_points[i + 1]
-        area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
-      next_point = horizontal_points[0]
-      point = horizontal_points[len(horizontal_points) - 1]
-      area += (next_point[1] + point[1]) / 2 * (next_point[0] - point[0])
-      self._area = abs(area)
+      self._area = 0
+      for triangle in self.triangles:
+        ab = np.zeros(3)
+        ac = np.zeros(3)
+        for i in range(0, 3):
+          ab[i] = triangle.coordinates[1][i] - triangle.coordinates[0][i]
+          ac[i] = triangle.coordinates[2][i] - triangle.coordinates[0][i]
+        self._area += np.linalg.norm(np.cross(ab, ac)) / 2
     return self._area
 
-  @property
-  def _points_rotated_to_horizontal(self):
-    """
-    polygon points rotated to horizontal
-    :return: [float]
-    """
-    z_vector = [0, 0, 1]
-    normal_vector = self.normal
-    horizontal_points = []
-    x = normal_vector[0]
-    y = normal_vector[1]
-
-    if x == 0 and y == 0:
-      # Already horizontal
-      for point in self.points:
-        horizontal_points.append([point.coordinates[0], point.coordinates[1], 0])
-    else:
-      alpha = self._angle_between_vectors(normal_vector, z_vector)
-      rotation_line = np.cross(normal_vector, z_vector)
-      third_axis = np.cross(normal_vector, rotation_line)
-      w_1 = rotation_line / np.linalg.norm(rotation_line)
-      w_2 = normal_vector
-      w_3 = third_axis / np.linalg.norm(third_axis)
-      rotation_matrix = np.array([[1, 0, 0],
-                                  [0, np.cos(alpha), -np.sin(alpha)],
-                                  [0, np.sin(alpha), np.cos(alpha)]])
-      base_matrix = np.array([w_1, w_2, w_3])
-      rotation_base_matrix = np.matmul(base_matrix.transpose(), rotation_matrix.transpose())
-      rotation_base_matrix = np.matmul(rotation_base_matrix, base_matrix)
-
-      if rotation_base_matrix is None:
-        sys.stderr.write('Warning: rotation base matrix returned None\n')
-      else:
-        for point in self.points:
-          new_point = np.matmul(rotation_base_matrix, point.coordinates)
-          horizontal_points.append(new_point)
-    return horizontal_points
-
   @property
   def normal(self) -> np.ndarray:
     """
@@ -275,283 +214,74 @@ class Polygon:
       return alpha
     return -alpha
 
-  def triangulate(self) -> List[Polygon]:
-    """
-    Triangulates a polygon following the ear clipping methodology
-    :return: list[triangles]
-    """
-    # todo: review triangulate_polygon in
-    #  https://github.com/mikedh/trimesh/blob/dad11126742e140ef46ba12f8cb8643c83356467/trimesh/creation.py#L415,
-    #  it had a problem with a class called 'triangle', but, if solved,
-    #  it could be a very good substitute of this method
-    #  this method is very dirty and has an infinite loop solved with a counter!!
-    if self._triangles is None:
-      points_list = self.points_list
-      normal = self.normal
-      if np.linalg.norm(normal) == 0:
-        sys.stderr.write('Not able to triangulate polygon\n')
-        return [self]
-      # are points concave or convex?
-      total_points_list, concave_points, convex_points = self._starting_lists(points_list, normal)
-
-      # list of ears
-      ears = []
-      j = 0
-      while (len(concave_points) > 3 or len(convex_points) != 0) and j < 100:
-        j += 1
-        for i in range(0, len(concave_points)):
-          ear = self._triangle(points_list, total_points_list, concave_points[i])
-          rest_points = []
-          for points in total_points_list:
-            rest_points.append(list(self.coordinates[points]))
-          if self._is_ear(ear, rest_points):
-            ears.append(ear)
-            point_to_remove = concave_points[i]
-            previous_point_in_list, next_point_in_list = self._enveloping_points(point_to_remove,
-                                                                                 total_points_list)
-            total_points_list.remove(point_to_remove)
-            concave_points.remove(point_to_remove)
-            # Was any of the adjacent points convex? -> check if changed status to concave
-            for convex_point in convex_points:
-              if convex_point == previous_point_in_list:
-                concave_points, convex_points, end_loop = self._if_concave_change_status(normal,
-                                                                                         points_list,
-                                                                                         convex_point,
-                                                                                         total_points_list,
-                                                                                         concave_points,
-                                                                                         convex_points,
-                                                                                         previous_point_in_list)
-                if end_loop:
-                  break
-                continue
-              if convex_point == next_point_in_list:
-                concave_points, convex_points, end_loop = self._if_concave_change_status(normal,
-                                                                                         points_list,
-                                                                                         convex_point,
-                                                                                         total_points_list,
-                                                                                         concave_points,
-                                                                                         convex_points,
-                                                                                         next_point_in_list)
-                if end_loop:
-                  break
-                continue
-            break
-        if len(total_points_list) <= 3 and len(convex_points) > 0:
-          sys.stderr.write('Not able to triangulate polygon\n')
-          return [self]
-      if j >= 100:
-        sys.stderr.write('Not able to triangulate polygon\n')
-        return [self]
-      last_ear = self._triangle(points_list, total_points_list, concave_points[1])
-      ears.append(last_ear)
-      self._triangles = ears
-    return self._triangles
-
   @staticmethod
-  def _starting_lists(points_list, normal) -> [List[float], List[float], List[float]]:
-    """
-    creates the list of vertices (points) that define the polygon (total_points_list), together with other two lists
-    separating points between convex and concave
-    :param points_list: points_list
-    :param normal: normal
-    :return: list[point], list[point], list[point]
-    """
-    concave_points = []
-    convex_points = []
-    # lists of concave and convex points
-    # case 1: first point
-    point = points_list[0:3]
-    previous_point = points_list[len(points_list) - 3:]
-    next_point = points_list[3:6]
-    index = 0
-    total_points_list = [index]
-    if Polygon._point_is_concave(normal, point, previous_point, next_point):
-      concave_points.append(index)
-    else:
-      convex_points.append(index)
-    # case 2: all points except first and last
-    for i in range(0, int((len(points_list) - 6) / 3)):
-      point = points_list[(i + 1) * 3:(i + 2) * 3]
-      previous_point = points_list[i * 3:(i + 1) * 3]
-      next_point = points_list[(i + 2) * 3:(i + 3) * 3]
-      index = i + 1
-      total_points_list.append(index)
-      if Polygon._point_is_concave(normal, point, previous_point, next_point):
-        concave_points.append(index)
-      else:
-        convex_points.append(index)
-    # case 3: last point
-    point = points_list[len(points_list) - 3:]
-    previous_point = points_list[len(points_list) - 6:len(points_list) - 3]
-    next_point = points_list[0:3]
-    index = int(len(points_list) / 3) - 1
-    total_points_list.append(index)
-    if Polygon._point_is_concave(normal, point, previous_point, next_point):
-      concave_points.append(index)
-    else:
-      convex_points.append(index)
-    return total_points_list, concave_points, convex_points
+  def triangle_mesh(vertices, normal):
+    min_x = 1e16
+    min_y = 1e16
+    min_z = 1e16
+    for vertex in vertices:
+      if vertex[0] < min_x:
+        min_x = vertex[0]
+      if vertex[1] < min_y:
+        min_y = vertex[1]
+      if vertex[2] < min_z:
+        min_z = vertex[2]
 
-  @staticmethod
-  def _triangle(points_list, total_points_list, point_position) -> Polygon:
-    """
-    creates a triangular polygon out of three points
-    :param points_list: points_list
-    :param total_points_list: [point]
-    :param point_position: int
-    :return: polygon
-    """
-    index = point_position * 3
-    previous_point_index, next_point_index = Polygon._enveloping_points_indices(point_position, total_points_list)
-    points = points_list[previous_point_index:previous_point_index + 3]
-    points = np.append(points, points_list[index:index + 3])
-    points = np.append(points, points_list[next_point_index:next_point_index + 3])
-    rows = points.size // 3
-    points = points.reshape(rows, 3)
-    triangle = Polygon(points)
-    return triangle
+    new_vertices = []
+    for vertex in vertices:
+      vertex = [vertex[0]-min_x, vertex[1]-min_y, vertex[2]-min_z]
+      new_vertices.append(vertex)
 
-  @staticmethod
-  def _enveloping_points_indices(point_position, total_points_list):
-    """
-    due to the fact that the lists are not circular, a method to find the previous and next points
-    of an specific one is needed
-    :param point_position: int
-    :param total_points_list: [point]
-    :return: int, int
-    """
-    previous_point_index = None
-    next_point_index = None
-    if point_position == total_points_list[0]:
-      previous_point_index = total_points_list[len(total_points_list) - 1] * 3
-      next_point_index = total_points_list[1] * 3
-    if point_position == total_points_list[len(total_points_list) - 1]:
-      previous_point_index = total_points_list[len(total_points_list) - 2] * 3
-      next_point_index = total_points_list[0] * 3
-    for i in range(1, len(total_points_list) - 1):
-      if point_position == total_points_list[i]:
-        previous_point_index = total_points_list[i - 1] * 3
-        next_point_index = total_points_list[i + 1] * 3
-    return previous_point_index, next_point_index
+    transformation_matrix = trimesh.geometry.plane_transform(origin=new_vertices[0], normal=normal)
 
-  @staticmethod
-  def _enveloping_points(point_to_remove, total_points_list):
-    """
-    due to the fact that the lists are not circular, a method to find the previous and next points
-    of an specific one is needed
-    :param point_to_remove: point
-    :param total_points_list: [point]
-    :return: point, point
-    """
-    index = total_points_list.index(point_to_remove)
-    if index == 0:
-      previous_point_in_list = total_points_list[len(total_points_list) - 1]
-      next_point_in_list = total_points_list[1]
-    elif index == len(total_points_list) - 1:
-      previous_point_in_list = total_points_list[len(total_points_list) - 2]
-      next_point_in_list = total_points_list[0]
-    else:
-      previous_point_in_list = total_points_list[index - 1]
-      next_point_in_list = total_points_list[index + 1]
-    return previous_point_in_list, next_point_in_list
+    coordinates = []
+    for vertex in vertices:
+      transformed_vertex = [vertex[0]-min_x, vertex[1]-min_y, vertex[2]-min_z, 1]
+      transformed_vertex = np.dot(transformation_matrix, transformed_vertex)
+      coordinate = [transformed_vertex[0], transformed_vertex[1]]
+      coordinates.append(coordinate)
 
-  @staticmethod
-  def _is_ear(ear, points) -> bool:
-    """
-    finds whether a triangle is an ear of the polygon
-    :param ear: polygon
-    :param points: [point]
-    :return: boolean
-    """
-    area_ear = ear.area
-    for point in points:
-      area_points = 0
-      point_is_not_vertex = True
+    polygon = shapley_polygon(coordinates)
+
+    try:
+      vertices_2d, faces = trimesh.creation.triangulate_polygon(polygon, engine='triangle')
+      mesh = Trimesh(vertices=vertices, faces=faces)
+
+      # check orientation
+      normal_sum = 0
       for i in range(0, 3):
-        if abs(np.linalg.norm(point) - np.linalg.norm(ear.coordinates[i])) < 0.0001:
-          point_is_not_vertex = False
-          break
-      if point_is_not_vertex:
-        for i in range(0, 3):
-          if i != 2:
-            new_points = ear.coordinates[i][:]
-            new_points = np.append(new_points, ear.coordinates[i + 1][:])
-            new_points = np.append(new_points, point[:])
-          else:
-            new_points = ear.coordinates[i][:]
-            new_points = np.append(new_points, point[:])
-            new_points = np.append(new_points, ear.coordinates[0][:])
-          rows = new_points.size // 3
-          new_points = new_points.reshape(rows, 3)
-          new_triangle = Polygon(new_points)
-          area_points += new_triangle.area
-      if abs(area_points - area_ear) < 1e-6:
-        # point_inside_ear = True
-        return False
-    return True
+        normal_sum += normal[i] + mesh.face_normals[0][i]
 
-  @staticmethod
-  def _if_concave_change_status(normal, points_list, convex_point, total_points_list,
-                                concave_points, convex_points, point_in_list) -> [List[float], List[float], bool]:
-    """
-    checks whether an convex specific point change its status to concave after removing one ear in the polygon
-    returning the new convex and concave points lists together with a flag advising that the list of total points
-    already 3 and, therefore, the triangulation must be finished.
-    :param normal: normal
-    :param points_list: points_list
-    :param convex_point: int
-    :param total_points_list: [point]
-    :param concave_points: [point]
-    :param convex_points: [point]
-    :param point_in_list: int
-    :return: list[points], list[points], boolean
-    """
-    end_loop = False
-    point = points_list[point_in_list * 3:(point_in_list + 1) * 3]
-    pointer = total_points_list.index(point_in_list) - 1
-    if pointer < 0:
-      pointer = len(total_points_list) - 1
-    previous_point = points_list[total_points_list[pointer] * 3:total_points_list[pointer] * 3 + 3]
-    pointer = total_points_list.index(point_in_list) + 1
-    if pointer >= len(total_points_list):
-      pointer = 0
-    next_point = points_list[total_points_list[pointer] * 3:total_points_list[pointer] * 3 + 3]
-    if Polygon._point_is_concave(normal, point, previous_point, next_point):
-      if concave_points[0] > convex_point:
-        concave_points.insert(0, convex_point)
-      elif concave_points[len(concave_points) - 1] < convex_point:
-        concave_points.append(convex_point)
-      else:
-        for point_index in range(0, len(concave_points) - 1):
-          if concave_points[point_index] < convex_point < concave_points[point_index + 1]:
-            concave_points.insert(point_index + 1, convex_point)
-      convex_points.remove(convex_point)
-      end_loop = True
-    return concave_points, convex_points, end_loop
+      if abs(normal_sum) <= 1E-10:
+        new_faces = []
+        for face in faces:
+          new_face = []
+          for i in range(0, len(face)):
+            new_face.append(face[len(face)-i-1])
+          new_faces.append(new_face)
+        mesh = Trimesh(vertices=vertices, faces=new_faces)
 
-  @staticmethod
-  def _point_is_concave(normal, point, previous_point, next_point) -> bool:
-    """
-    returns whether a point is concave
-    :param normal: normal
-    :param point: point
-    :param previous_point: point
-    :param next_point: point
-    :return: boolean
-    """
-    is_concave = False
-    accepted_error = 0.1
-    points = np.append(previous_point, point)
-    points = np.append(points, next_point)
-    rows = points.size // 3
-    points = points.reshape(rows, 3)
-    triangle = Polygon(points)
-    error_sum = 0
-    for i in range(0, len(normal)):
-      error_sum += triangle.normal[i] - normal[i]
-    if np.abs(error_sum) < accepted_error:
-      is_concave = True
-    return is_concave
+      return mesh
+
+    except ValueError:
+      logger.error(f'Not able to triangulate polygon\n')
+      sys.stderr.write(f'Not able to triangulate polygon\n')
+      _vertices = [[0, 0, 0], [0, 0, 1], [0, 1, 0]]
+      _faces = [[0, 1, 2]]
+      return Trimesh(vertices=_vertices, faces=_faces)
+
+  @property
+  def triangles(self) -> List[Polygon]:
+    if self._triangles is None:
+      self._triangles = []
+      _mesh = self.triangle_mesh(self.coordinates, self.normal)
+      for face in _mesh.faces:
+        points = []
+        for vertex in face:
+          points.append(self.coordinates[vertex])
+        polygon = Polygon(points)
+        self._triangles.append(polygon)
+    return self._triangles
 
   @staticmethod
   def _angle_between_vectors(vec_1, vec_2):
@@ -652,12 +382,12 @@ class Polygon:
   @property
   def vertices(self) -> np.ndarray:
     """
-    Get polyhedron vertices
+    Get polygon vertices
     :return: np.ndarray(int)
     """
     if self._vertices is None:
       vertices, self._vertices = [], []
-      _ = [vertices.extend(s.coordinates) for s in self.triangulate()]
+      _ = [vertices.extend(s.coordinates) for s in self.triangles]
       for vertex_1 in vertices:
         found = False
         for vertex_2 in self._vertices:
@@ -677,17 +407,17 @@ class Polygon:
   @property
   def faces(self) -> List[List[int]]:
     """
-    Get polyhedron triangular faces
+    Get polygon triangular faces
     :return: [face]
     """
     if self._faces is None:
       self._faces = []
 
-      for polygon in self.triangulate():
+      for polygon in self.triangles:
         face = []
         points = polygon.coordinates
         if len(points) != 3:
-          sub_polygons = polygon.triangulate()
+          sub_polygons = polygon.triangles
           # todo: I modified this! To be checked @Guille
           if len(sub_polygons) >= 1:
             for sub_polygon in sub_polygons:
diff --git a/hub/city_model_structure/attributes/polyhedron.py b/hub/city_model_structure/attributes/polyhedron.py
index f27a8ff1..553f1e26 100644
--- a/hub/city_model_structure/attributes/polyhedron.py
+++ b/hub/city_model_structure/attributes/polyhedron.py
@@ -91,7 +91,7 @@ class Polyhedron:
         face = []
         points = polygon.coordinates
         if len(points) != 3:
-          sub_polygons = polygon.triangulate()
+          sub_polygons = polygon.triangles
           # todo: I modified this! To be checked @Guille
           if len(sub_polygons) >= 1:
             for sub_polygon in sub_polygons:
diff --git a/hub/city_model_structure/building.py b/hub/city_model_structure/building.py
index 296182f1..d5129142 100644
--- a/hub/city_model_structure/building.py
+++ b/hub/city_model_structure/building.py
@@ -6,8 +6,10 @@ Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
 Code contributors: Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
 """
 
+import sys
 from typing import List, Union
 import numpy as np
+from hub.hub_logger import logger
 import hub.helpers.constants as cte
 from hub.city_model_structure.building_demand.surface import Surface
 from hub.city_model_structure.city_object import CityObject
@@ -46,20 +48,31 @@ class Building(CityObject):
     self._roofs = []
     self._walls = []
     self._internal_walls = []
+    self._ground_walls = []
+    self._attic_floors = []
+    self._interior_slabs = []
     for surface_id, surface in enumerate(self.surfaces):
       self._min_x = min(self._min_x, surface.lower_corner[0])
       self._min_y = min(self._min_y, surface.lower_corner[1])
       self._min_z = min(self._min_z, surface.lower_corner[2])
       surface.id = surface_id
-      # todo: consider all type of surfaces, not only these four
       if surface.type == cte.GROUND:
         self._grounds.append(surface)
       elif surface.type == cte.WALL:
         self._walls.append(surface)
       elif surface.type == cte.ROOF:
         self._roofs.append(surface)
-      else:
+      elif surface.type == cte.INTERIOR_WALL:
         self._internal_walls.append(surface)
+      elif surface.type == cte.GROUND_WALL:
+        self._ground_walls.append(surface)
+      elif surface.type == cte.ATTIC_FLOOR:
+        self._attic_floors.append(surface)
+      elif surface.type == cte.INTERIOR_SLAB:
+        self._interior_slabs.append(surface)
+      else:
+        logger.error(f'Building {self.name} [alias {self.alias}]  has an unexpected surface type {surface.type}.\n')
+        sys.stderr.write(f'Building {self.name} [alias {self.alias}]  has an unexpected surface type {surface.type}.\n')
 
   @property
   def shell(self) -> Polyhedron:
diff --git a/hub/city_model_structure/building_demand/occupancy.py b/hub/city_model_structure/building_demand/occupancy.py
index 99e4ab23..46d34aa1 100644
--- a/hub/city_model_structure/building_demand/occupancy.py
+++ b/hub/city_model_structure/building_demand/occupancy.py
@@ -6,7 +6,6 @@ Project Coder Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
 """
 from typing import Union, List
 from hub.city_model_structure.attributes.schedule import Schedule
-from hub.city_model_structure.building_demand.occupant import Occupant
 
 
 class Occupancy:
@@ -106,19 +105,3 @@ class Occupancy:
     :param value: [Schedule]
     """
     self._occupancy_schedules = value
-
-  @property
-  def occupants(self) -> Union[None, List[Occupant]]:
-    """
-    Get list of occupants
-    :return: None or List of Occupant
-    """
-    return self._occupants
-
-  @occupants.setter
-  def occupants(self, value):
-    """
-    Set list of occupants
-    :param value: [Occupant]
-    """
-    self._occupants = value
diff --git a/hub/city_model_structure/building_demand/occupant.py b/hub/city_model_structure/building_demand/occupant.py
deleted file mode 100644
index a5916d80..00000000
--- a/hub/city_model_structure/building_demand/occupant.py
+++ /dev/null
@@ -1,145 +0,0 @@
-"""
-Occupant module
-SPDX - License - Identifier: LGPL - 3.0 - or -later
-Copyright © 2022 Concordia CERC group
-Project Coder Sanam Dabirian sanam.dabirian@mail.concordia.ca
-Code contributors: Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
-"""
-
-
-class Occupant:
-  """
-  Occupant class
-  """
-
-  def __init__(self):
-    """
-    Constructor
-    """
-
-    self._heat_dissipation = None
-    self._occupancy_rate = None
-    self._occupant_type = None
-    self._arrival_time = None
-    self._departure_time = None
-    self._break_time = None
-    self._day_of_week = None
-    self._pd_of_meetings_duration = None
-
-  @property
-  def heat_dissipation(self):
-    """
-    Get heat dissipation of occupants in W/person
-    :return: float
-    """
-    return self._heat_dissipation
-
-  @heat_dissipation.setter
-  def heat_dissipation(self, value):
-    """
-    Set heat dissipation of occupants in W/person
-    :param value: float
-    """
-    self._heat_dissipation = float(value)
-
-  @property
-  def occupancy_rate(self):
-    """
-    Get rate of schedules
-    :return: float
-    """
-    return self._occupancy_rate
-
-  @occupancy_rate.setter
-  def occupancy_rate(self, value):
-    """
-    Set rate of schedules
-    :param value: float
-    """
-    self._occupancy_rate = float(value)
-
-  @property
-  def occupant_type(self):
-    """
-    Get type of schedules
-    :return: str
-    """
-    return self._occupant_type
-
-  @occupant_type.setter
-  def occupant_type(self, value):
-    """
-    Set type of schedules
-    :param value: float
-    """
-    self._occupant_type = float(value)
-
-  @property
-  def arrival_time(self):
-    """
-    Get the arrival time of the occupant (for office building) in UTC with format YYYYMMDD HH:mm:ss
-    :return: time
-    """
-    return self._arrival_time
-
-  @arrival_time.setter
-  def arrival_time(self, value):
-    """
-    Set the arrival time of the occupant (for office building) in UTC with format YYYYMMDD HH:mm:ss
-    :param value: time
-    """
-    self._arrival_time = value
-
-  @property
-  def departure_time(self):
-    """
-    Get the departure time of the occupant (for office building) in UTC with format YYYYMMDD HH:mm:ss
-    :return: time
-    """
-    return self._departure_time
-
-  @departure_time.setter
-  def departure_time(self, value):
-    """
-    Set the departure time of the occupant (for office building) in UTC with format YYYYMMDD HH:mm:ss
-    :param value: str
-    """
-    self._departure_time = value
-
-  @property
-  def break_time(self):
-    """
-    Get the lunch or break time of the occupant (for office building) in UTC with format ????
-    :return: break time
-    """
-    # todo @Sanam: define this format, is it the starting time? is it a list with both, starting and ending time?
-    return self._break_time
-
-  @property
-  def day_of_week(self):
-    """
-    Get the day of the week (MON, TUE, WED, THU, FRI, SAT, SUN)
-    :return: str
-    """
-    # todo @Sanam: is this a property or should it be a function
-    #  to get the day of the week of an specific day of the year?
-    return self._day_of_week
-
-  @property
-  def pd_of_meetings_duration(self):
-    """
-    Get the probability distribution of the meeting duration
-    :return: ??
-    """
-    # todo @Sanam: what format are you expecting here??
-    return self._pd_of_meetings_duration
-
-  @pd_of_meetings_duration.setter
-  def pd_of_meetings_duration(self, value):
-    """
-    Get the probability distribution of the meeting duration
-    :param value: ??
-    :return:
-    """
-    # todo @Sanam: what format are you expecting here??
-    self._pd_of_meetings_duration = value
diff --git a/hub/city_model_structure/building_demand/surface.py b/hub/city_model_structure/building_demand/surface.py
index 1c4e4cac..dc392160 100644
--- a/hub/city_model_structure/building_demand/surface.py
+++ b/hub/city_model_structure/building_demand/surface.py
@@ -72,14 +72,6 @@ class Surface:
     if value is not None:
       self._id = str(value)
 
-  # todo: implement share surfaces
-  @property
-  def share_surfaces(self):
-    """
-    Raises not implemented error
-    """
-    raise NotImplementedError
-
   def _max_coord(self, axis):
     if axis == 'x':
       axis = 0
@@ -163,11 +155,11 @@ class Surface:
   @property
   def type(self):
     """
-    Get surface type Ground, Wall or Roof
+    Get surface type Ground, Ground wall, Wall, Attic floor, Interior slab, Interior wall, Roof or Virtual internal
+    If the geometrical LoD is lower than 4,
+    the surfaces' types are not defined in the importer and can only be Ground, Wall or Roof
     :return: str
     """
-
-    # todo: there are more types: internal wall, internal floor... this method must be redefined
     if self._type is None:
       grad = np.rad2deg(self.inclination)
       if grad >= 170:
@@ -304,7 +296,7 @@ class Surface:
     :return: Surface, Surface, Any
     """
     # todo: check return types
-    # todo: recheck this method for LoD3 (windows)
+    #  recheck this method for LoD3 (windows)
     origin = Point([0, 0, z])
     normal = np.array([0, 0, 1])
     plane = Plane(normal=normal, origin=origin)
diff --git a/hub/city_model_structure/city.py b/hub/city_model_structure/city.py
index ef8ebbd1..be5d954f 100644
--- a/hub/city_model_structure/city.py
+++ b/hub/city_model_structure/city.py
@@ -23,7 +23,6 @@ from hub.city_model_structure.iot.station import Station
 from hub.city_model_structure.level_of_detail import LevelOfDetail
 from hub.city_model_structure.machine import Machine
 from hub.city_model_structure.parts_consisting_building import PartsConsistingBuilding
-from hub.city_model_structure.subway_entrance import SubwayEntrance
 from hub.helpers.geometry_helper import GeometryHelper
 from hub.helpers.location import Location
 from hub.city_model_structure.energy_system import EnergySystem
@@ -40,9 +39,7 @@ class City:
     self._lower_corner = lower_corner
     self._upper_corner = upper_corner
     self._buildings = None
-    self._subway_entrances = None
     self._srs_name = srs_name
-    # todo: right now extracted at city level, in the future should be extracted also at building level if exist
     self._location = None
     self._country_code = None
     self._climate_reference_city = None
@@ -163,9 +160,6 @@ class City:
       if self.buildings is not None:
         for building in self.buildings:
           self._city_objects.append(building)
-      if self.subway_entrances is not None:
-        for subway_entrance in self.subway_entrances:
-          self._city_objects.append(subway_entrance)
       if self.energy_systems is not None:
         for energy_system in self.energy_systems:
           self._city_objects.append(energy_system)
@@ -179,14 +173,6 @@ class City:
     """
     return self._buildings
 
-  @property
-  def subway_entrances(self) -> Union[List[SubwayEntrance], None]:
-    """
-    Get the subway entrances belonging to the city
-    :return: a list of subway entrances objects or none
-    """
-    return self._subway_entrances
-
   @property
   def lower_corner(self) -> List[float]:
     """
@@ -224,10 +210,6 @@ class City:
       if self._buildings is None:
         self._buildings = []
       self._buildings.append(new_city_object)
-    elif new_city_object.type == 'subway_entrance':
-      if self._subway_entrances is None:
-        self._subway_entrances = []
-      self._subway_entrances.append(new_city_object)
     elif new_city_object.type == 'energy_system':
       if self._energy_systems is None:
         self._energy_systems = []
diff --git a/hub/city_model_structure/subway_entrance.py b/hub/city_model_structure/subway_entrance.py
deleted file mode 100644
index cd0c8598..00000000
--- a/hub/city_model_structure/subway_entrance.py
+++ /dev/null
@@ -1,45 +0,0 @@
-"""
-Subway entrance module
-SPDX - License - Identifier: LGPL - 3.0 - or -later
-Copyright © 2022 Concordia CERC group
-Project Coder Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
-"""
-from hub.city_model_structure.city_object import CityObject
-
-
-class SubwayEntrance(CityObject):
-  """
-  SubwayEntrance(CityObject) class
-  """
-  def __init__(self, name, latitude, longitude):
-    super().__init__(name, 0)
-    self._name = name
-    self._latitude = latitude
-    self._longitude = longitude
-    self._type = 'subway_entrance'
-
-  @property
-  def latitude(self):
-    # todo: to be defined the spacial point and the units
-    """
-    Get latitude
-    :return: float
-    """
-    return self._latitude
-
-  @property
-  def longitude(self):
-    # todo: to be defined the spacial point and the units
-    """
-    Get longitude
-    :return: float
-    """
-    return self._longitude
-
-  @property
-  def name(self):
-    """
-    Get name
-    :return: str
-    """
-    return self._name
diff --git a/hub/helpers/geometry_helper.py b/hub/helpers/geometry_helper.py
index 1c67bf4a..24c202e3 100644
--- a/hub/helpers/geometry_helper.py
+++ b/hub/helpers/geometry_helper.py
@@ -39,59 +39,12 @@ class GeometryHelper:
     max_distance = ConfigurationHelper().max_location_distance_for_shared_walls
     return GeometryHelper.distance_between_points(location1, location2) < max_distance
 
-  def almost_same_area(self, area_1, area_2):
-    """
-    Compare two areas and decides if they are almost equal (absolute error under delta)
-    :param area_1
-    :param area_2
-    :return: Boolean
-    """
-    if area_1 == 0 or area_2 == 0:
-      return False
-    delta = math.fabs(area_1 - area_2)
-    return delta <= self._area_delta
-
-  def is_almost_same_surface(self, surface_1, surface_2):
-    """
-    Compare two surfaces and decides if they are almost equal (quadratic error under delta)
-    :param surface_1: Surface
-    :param surface_2: Surface
-    :return: Boolean
-    """
-
-    # delta is grads an need to be converted into radians
-    delta = np.rad2deg(self._delta)
-    difference = (surface_1.inclination - surface_2.inclination) % math.pi
-    if abs(difference) > delta:
-      return False
-    # s1 and s2 are at least almost parallel surfaces
-    # calculate distance point to plane using all the vertex
-    # select surface1 value for the point (X,Y,Z) where two of the values are 0
-    minimum_distance = self._delta + 1
-    parametric = surface_2.polygon.get_parametric()
-    normal_2 = surface_2.normal
-    for point in surface_1.points:
-      distance = abs(
-        (point[0] * parametric[0]) + (point[1] * parametric[1]) + (point[2] * parametric[2]) + parametric[3])
-      normal_module = math.sqrt(pow(normal_2[0], 2) + pow(normal_2[1], 2) + pow(normal_2[2], 2))
-
-      if normal_module == 0:
-        continue
-      distance = distance / normal_module
-      if distance < minimum_distance:
-        minimum_distance = distance
-      if minimum_distance <= self._delta:
-        break
-
-    if minimum_distance > self._delta or surface_1.intersect(surface_2) is None:
-      return False
-    return True
-
   @staticmethod
   def segment_list_to_trimesh(lines) -> Trimesh:
     """
     Transform a list of segments into a Trimesh
     """
+    # todo: trimesh has a method for this
     line_points = [lines[0][0], lines[0][1]]
     lines.remove(lines[0])
     while len(lines) > 1:
@@ -106,7 +59,7 @@ class GeometryHelper:
           line_points.append(line[0])
           lines.pop(i - 1)
           break
-    polyhedron = Polyhedron(Polygon(line_points).triangulate())
+    polyhedron = Polyhedron(Polygon(line_points).triangles)
     trimesh = Trimesh(polyhedron.vertices, polyhedron.faces)
     return trimesh
 
diff --git a/hub/helpers/monthly_to_hourly_demand.py b/hub/helpers/monthly_to_hourly_demand.py
deleted file mode 100644
index a83b9673..00000000
--- a/hub/helpers/monthly_to_hourly_demand.py
+++ /dev/null
@@ -1,138 +0,0 @@
-"""
-monthly_to_hourly_demand module
-SPDX - License - Identifier: LGPL - 3.0 - or -later
-Copyright © 2022 Concordia CERC group
-Project Coder Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
-"""
-import calendar as cal
-import pandas as pd
-from hub.city_model_structure.building_demand.occupant import Occupant
-import hub.helpers.constants as cte
-
-
-class MonthlyToHourlyDemand:
-  """
-  MonthlyToHourlyDemand class
-  """
-  def __init__(self, building, conditioning_seasons):
-    self._hourly_heating = pd.DataFrame()
-    self._hourly_cooling = pd.DataFrame()
-    self._building = building
-    self._conditioning_seasons = conditioning_seasons
-
-  def hourly_heating(self, key):
-    """
-    hourly distribution of the monthly heating of a building
-    :param key: string
-    :return: [hourly_heating]
-    """
-    # todo: this method and the insel model have to be reviewed for more than one thermal zone
-    external_temp = self._building.external_temperature[cte.HOUR]
-    # todo: review index depending on how the schedules are defined, either 8760 or 24 hours
-    for usage in self._building.usages:
-      temp_set = float(usage.heating_setpoint)-3
-      temp_back = float(usage.heating_setback)-3
-      # todo: if these are data frames, then they should be called as (Occupancy should be in low case):
-      # usage.schedules.Occupancy
-      # self._conditioning_seasons.heating
-      occupancy = Occupant().get_complete_year_schedule(usage.schedules['Occupancy'])
-      heating_schedule = self._conditioning_seasons['heating']
-
-      hourly_heating = []
-      i = 0
-      j = 0
-      temp_grad_day = []
-      for month in range(1, 13):
-        temp_grad_month = 0
-        month_range = cal.monthrange(2015, month)[1]
-        for _ in range(1, month_range+1):
-          external_temp_med = 0
-          for hour in range(0, 24):
-            external_temp_med += external_temp[key][i]/24
-          for hour in range(0, 24):
-            if external_temp_med < temp_set and heating_schedule[month-1] == 1:
-              if occupancy[hour] > 0:
-                hdd = temp_set - external_temp[key][i]
-                if hdd < 0:
-                  hdd = 0
-                temp_grad_day.append(hdd)
-              else:
-                hdd = temp_back - external_temp[key][i]
-                if hdd < 0:
-                  hdd = 0
-                temp_grad_day.append(hdd)
-            else:
-              temp_grad_day.append(0)
-
-            temp_grad_month += temp_grad_day[i]
-            i += 1
-
-        for _ in range(1, month_range + 1):
-          for hour in range(0, 24):
-            monthly_demand = self._building.heating[cte.MONTH][month-1]
-            if monthly_demand == 'NaN':
-              monthly_demand = 0
-            if temp_grad_month == 0:
-              hourly_demand = 0
-            else:
-              hourly_demand = float(monthly_demand)*float(temp_grad_day[j])/float(temp_grad_month)
-            hourly_heating.append(hourly_demand)
-            j += 1
-      self._hourly_heating = pd.DataFrame(data=hourly_heating, columns=['monthly to hourly'])
-    return self._hourly_heating
-
-  def hourly_cooling(self, key):
-    """
-    hourly distribution of the monthly cooling of a building
-    :param key: string
-    :return: [hourly_cooling]
-    """
-    # todo: this method and the insel model have to be reviewed for more than one thermal zone
-    external_temp = self._building.external_temperature[cte.HOUR]
-    # todo: review index depending on how the schedules are defined, either 8760 or 24 hours
-    for usage in self._building.usages:
-      temp_set = float(usage.cooling_setpoint)
-      temp_back = 100
-      occupancy = Occupant().get_complete_year_schedule(usage.schedules['Occupancy'])
-      cooling_schedule = self._conditioning_seasons['cooling']
-
-      hourly_cooling = []
-      i = 0
-      j = 0
-      temp_grad_day = []
-      for month in range(1, 13):
-        temp_grad_month = 0
-        month_range = cal.monthrange(2015, month)[1]
-        for _ in range(1, month_range[1] + 1):
-          for hour in range(0, 24):
-            if external_temp[key][i] > temp_set and cooling_schedule[month - 1] == 1:
-              if occupancy[hour] > 0:
-                cdd = external_temp[key][i] - temp_set
-                if cdd < 0:
-                  cdd = 0
-                temp_grad_day.append(cdd)
-              else:
-                cdd = external_temp[key][i] - temp_back
-                if cdd < 0:
-                  cdd = 0
-                temp_grad_day.append(cdd)
-            else:
-              temp_grad_day.append(0)
-
-            temp_grad_month += temp_grad_day[i]
-            i += 1
-
-        for _ in range(1, month_range[1] + 1):
-          for hour in range(0, 24):
-            #            monthly_demand = self._building.heating[cte.MONTH]['INSEL'][month-1]
-            monthly_demand = self._building.cooling[cte.MONTH][month - 1]
-            if monthly_demand == 'NaN':
-              monthly_demand = 0
-            if temp_grad_month == 0:
-              hourly_demand = 0
-            else:
-              hourly_demand = float(monthly_demand) * float(temp_grad_day[j]) / float(temp_grad_month)
-            hourly_cooling.append(hourly_demand)
-            j += 1
-      self._hourly_cooling = pd.DataFrame(data=hourly_cooling, columns=['monthly to hourly'])
-    return self._hourly_cooling
diff --git a/hub/imports/geometry/osm_subway.py b/hub/imports/geometry/osm_subway.py
deleted file mode 100644
index 8e0bd27b..00000000
--- a/hub/imports/geometry/osm_subway.py
+++ /dev/null
@@ -1,58 +0,0 @@
-"""
-OsmSubway module parses osm files and import the metro location into the city model structure
-SPDX - License - Identifier: LGPL - 3.0 - or -later
-Copyright © 2022 Concordia CERC group
-Project Coder Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
-"""
-
-import sys
-import xmltodict
-from pyproj import Transformer
-from hub.city_model_structure.city import City
-from hub.city_model_structure.subway_entrance import SubwayEntrance
-
-
-class OsmSubway:
-  """
-  Open street map subway
-  """
-  def __init__(self, path):
-    self._city = None
-    self._subway_entrances = []
-    with open(path) as osm:
-      self._osm = xmltodict.parse(osm.read(), force_list='tag')
-      for node in self._osm['osm']['node']:
-        if 'tag' not in node:
-          continue
-        for tag in node['tag']:
-          if '@v' not in tag:
-            continue
-          if tag['@v'] == 'subway_entrance':
-            subway_entrance = SubwayEntrance(node['@id'], node['@lat'], node['@lon'])
-            self._subway_entrances.append(subway_entrance)
-
-  @property
-  def city(self) -> City:
-    """
-    Get a city with subway entrances
-    """
-    transformer = Transformer.from_crs("EPSG:4326", "EPSG:3857")
-    lower_corner = [sys.float_info.max, sys.float_info.max, 0]
-    upper_corner = [sys.float_info.min, sys.float_info.min, 0]
-    x = 0
-    y = 1
-    for subway_entrance in self._subway_entrances:
-      coordinate = transformer.transform(subway_entrance.longitude, subway_entrance.latitude)
-      if coordinate[x] >= upper_corner[x]:
-        upper_corner[x] = coordinate[x]
-      if coordinate[y] >= upper_corner[y]:
-        upper_corner[y] = coordinate[y]
-      if coordinate[x] < lower_corner[x]:
-        lower_corner[x] = coordinate[x]
-      if coordinate[y] < lower_corner[y]:
-        lower_corner[y] = coordinate[y]
-
-    city = City(lower_corner, upper_corner, 'unknown')
-    for subway_entrance in self._subway_entrances:
-      city.add_city_object(subway_entrance)
-    return city
diff --git a/hub/imports/geometry_factory.py b/hub/imports/geometry_factory.py
index b8ba7e68..2c870085 100644
--- a/hub/imports/geometry_factory.py
+++ b/hub/imports/geometry_factory.py
@@ -9,7 +9,6 @@ import geopandas
 from hub.city_model_structure.city import City
 from hub.imports.geometry.citygml import CityGml
 from hub.imports.geometry.obj import Obj
-from hub.imports.geometry.osm_subway import OsmSubway
 from hub.imports.geometry.rhino import Rhino
 from hub.imports.geometry.gpandas import GPandas
 from hub.imports.geometry.geojson import Geojson
@@ -83,14 +82,6 @@ class GeometryFactory:
                    self._function_field,
                    self._function_to_hub).city
 
-  @property
-  def _osm_subway(self) -> City:
-    """
-    Enrich the city by using OpenStreetMap information as data source
-    :return: City
-    """
-    return OsmSubway(self._path).city
-
   @property
   def _rhino(self) -> City:
     """