diff --git a/city_model_structure/attributes/plane.py b/city_model_structure/attributes/plane.py
new file mode 100644
index 00000000..f653d6c7
--- /dev/null
+++ b/city_model_structure/attributes/plane.py
@@ -0,0 +1,57 @@
+"""
+Plane module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2020 Project Author Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+"""
+
+from typing import TypeVar
+import numpy as np
+
+Point = TypeVar('Point')
+
+
+class Plane:
+  """
+  Plane class
+  """
+
+  def __init__(self, origin=None, normal=None):
+    # todo: other options to define the plane:
+    #  by two lines
+    #  by three points
+    self._origin = origin
+    self._normal = normal
+    self._opposite_normal = None
+
+  @property
+  def origin(self) -> Point:
+    """
+    Point origin of the plane
+    return Point
+    """
+    if self._origin is None:
+      raise NotImplementedError
+    return self._origin
+
+  @property
+  def normal(self):
+    """
+    Plane normal [x, y, z]
+    return np.ndarray
+    """
+    if self._normal is None:
+      raise NotImplementedError
+    return self._normal
+
+  @property
+  def opposite_normal(self):
+    """
+    Plane normal in the opposite direction [x, y, z]
+    return np.ndarray
+    """
+    if self._opposite_normal is None:
+      coordinates = []
+      for coordinate in self.normal:
+        coordinates.append(-coordinate)
+      self._opposite_normal = np.array(coordinates)
+    return self._opposite_normal
diff --git a/city_model_structure/attributes/polygon.py b/city_model_structure/attributes/polygon.py
index c9082e75..7705a575 100644
--- a/city_model_structure/attributes/polygon.py
+++ b/city_model_structure/attributes/polygon.py
@@ -10,6 +10,8 @@ import sys
 import numpy as np
 import math
 from city_model_structure.attributes.point import Point
+from trimesh import Trimesh
+import trimesh.intersections
 
 
 class Polygon:
@@ -29,6 +31,9 @@ class Polygon:
     self._edges = None
 #    self._coordinates = None
     self._coordinates = coordinates
+    self._triangles = None
+    self._vertices = None
+    self._faces = None
 
 #  @property
 #  def points(self) -> List[Point]:
@@ -66,7 +71,7 @@ class Polygon:
   def coordinates(self) -> List[np.ndarray]:
     """
     List of points in the shape of its coordinates belonging to the polygon [[x, y, z],...]
-    :return: np.array
+    :return: [np.ndarray]
     """
     return self._coordinates
 
@@ -234,58 +239,60 @@ class Polygon:
     #  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!!
-    points_list = self.points_list
-    normal = self.normal
-    if np.linalg.norm(normal) == 0:
-      sys.stderr.write(f'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 p in total_points_list:
-          rest_points.append(list(self.coordinates[p]))
-        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:
+    if self._triangles is None:
+      points_list = self.points_list
+      normal = self.normal
+      if np.linalg.norm(normal) == 0:
         sys.stderr.write(f'Not able to triangulate polygon\n')
         return [self]
-    if j >= 100:
-      sys.stderr.write(f'Not able to triangulate polygon\n')
-      return [self]
-    last_ear = self._triangle(points_list, total_points_list, concave_points[1])
-    ears.append(last_ear)
-    return ears
+      # 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 p in total_points_list:
+            rest_points.append(list(self.coordinates[p]))
+          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(f'Not able to triangulate polygon\n')
+          return [self]
+      if j >= 100:
+        sys.stderr.write(f'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]]:
@@ -514,30 +521,42 @@ class Polygon:
 
   @property
   def inverse(self):
+    """
+    Flips the order of the coordinates
+    :return: [np.ndarray]
+    """
     if self._inverse is None:
-      self._inverse = self.points[::-1]
+      self._inverse = self.coordinates[::-1]
     return self._inverse
 
-#  def divide(self, polygon):
-    
-#    return polygon_1, polygon_2, intersection
+  def divide(self, plane):
+    """
+    Divides the polygon in two by a plane
+    :param plane: plane that intersects with self to divide it in two parts (Plane)
+    :return: Polygon, Polygon, [Point]
+    """
+    tri_polygons = Trimesh(vertices=self.vertices, faces=self.faces)
+    intersection = trimesh.intersections.mesh_plane(tri_polygons, plane.normal, plane.origin.coordinates)
+    polys_1 = trimesh.intersections.slice_mesh_plane(tri_polygons, plane.opposite_normal, plane.origin.coordinates)
+    polys_2 = trimesh.intersections.slice_mesh_plane(tri_polygons, plane.normal, plane.origin.coordinates)
+    triangles_1 = []
+    for triangle in polys_1.triangles:
+      triangles_1.append(Polygon(triangle))
+    polygon_1 = self._reshape(triangles_1)
+    triangles_2 = []
+    for triangle in polys_2.triangles:
+      triangles_2.append(Polygon(triangle))
+    polygon_2 = self._reshape(triangles_2)
+    return polygon_1, polygon_2, intersection
 
-  def reshape(self, triangles) -> Polygon:
+  def _reshape(self, triangles) -> Polygon:
     edges_list = []
     for i in range(0, len(triangles)):
       for edge in triangles[i].edges:
-        print('edge')
-        print(edge[0].coordinates, edge[1].coordinates)
         if not self._edge_in_edges_list(edge, edges_list):
           edges_list.append(edge)
-          print('list')
-          for e in edges_list:
-            print(e[0].coordinates, e[1].coordinates)
         else:
-          print('remove')
           edges_list = self._remove_from_list(edge, edges_list)
-          for e in edges_list:
-            print(e[0].coordinates, e[1].coordinates)
     points = self._order_points(edges_list)
     return Polygon(points)
 
@@ -551,16 +570,17 @@ class Polygon:
 
   @staticmethod
   def _order_points(edges_list):
+    # todo: not sure that this method works for any case -> RECHECK
     points = edges_list[0]
-    for i in range(1, len(edges_list)):
-      point_1 = edges_list[i][0]
-      point_2 = points[len(points)-1]
-      if point_1.distance_to_point(point_2) == 0:
-        points.append(edges_list[i][1])
+    for j in range(0, len(points)):
+      for i in range(1, len(edges_list)):
+        point_1 = edges_list[i][0]
+        point_2 = points[len(points)-1]
+        if point_1.distance_to_point(point_2) == 0:
+          points.append(edges_list[i][1])
     points.remove(points[len(points)-1])
     array_points = []
     for point in points:
-      print(point.coordinates)
       array_points.append(point.coordinates)
     return np.array(array_points)
 
@@ -572,3 +592,73 @@ class Polygon:
              (ed[1].distance_to_point(edge[0]) == 0 and ed[0].distance_to_point(edge[1]) == 0)):
         new_list.append(ed)
     return new_list
+
+  @property
+  def vertices(self) -> np.ndarray:
+    """
+    Polyhedron vertices
+    :return: np.ndarray(int)
+    """
+    if self._vertices is None:
+      vertices, self._vertices = [], []
+      _ = [vertices.extend(s.coordinates) for s in self.triangulate()]
+      for vertex_1 in vertices:
+        found = False
+        for vertex_2 in self._vertices:
+          found = False
+          power = 0
+          for dimension in range(0, 3):
+            power += math.pow(vertex_2[dimension] - vertex_1[dimension], 2)
+          distance = math.sqrt(power)
+          if distance == 0:
+            found = True
+            break
+        if not found:
+          self._vertices.append(vertex_1)
+      self._vertices = np.asarray(self._vertices)
+    return self._vertices
+
+  @property
+  def faces(self) -> List[List[int]]:
+    """
+    Polyhedron triangular faces
+    :return: [face]
+    """
+    if self._faces is None:
+      self._faces = []
+
+      for polygon in self.triangulate():
+        face = []
+        points = polygon.coordinates
+        if len(points) != 3:
+          sub_polygons = polygon.triangulate()
+          # todo: I modified this! To be checked @Guille
+          if len(sub_polygons) >= 1:
+            for sub_polygon in sub_polygons:
+              face = []
+              points = sub_polygon.coordinates
+              for point in points:
+                face.append(self._position_of(point, face))
+              self._faces.append(face)
+        else:
+          for point in points:
+            face.append(self._position_of(point, face))
+          self._faces.append(face)
+    return self._faces
+
+  def _position_of(self, point, face):
+    """
+    position of a specific point in the list of points that define a face
+    :return: int
+    """
+    vertices = self.vertices
+    for i in range(len(vertices)):
+      # ensure not duplicated vertex
+      power = 0
+      vertex2 = vertices[i]
+      for dimension in range(0, 3):
+        power += math.pow(vertex2[dimension] - point[dimension], 2)
+      distance = math.sqrt(power)
+      if i not in face and distance == 0:
+        return i
+    return -1
diff --git a/city_model_structure/attributes/storey.py b/city_model_structure/attributes/storey.py
new file mode 100644
index 00000000..568c506b
--- /dev/null
+++ b/city_model_structure/attributes/storey.py
@@ -0,0 +1,98 @@
+"""
+Storey module
+SPDX - License - Identifier: LGPL - 3.0 - or -later
+Copyright © 2020 Project Author Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
+"""
+
+from __future__ import annotations
+from typing import List
+import numpy as np
+
+from city_model_structure.attributes.surface import Surface
+from city_model_structure.attributes.thermal_boundary import ThermalBoundary
+from city_model_structure.attributes.thermal_zone import ThermalZone
+import helpers.constants as cte
+
+
+class Storey:
+  # todo: rethink this class for buildings with windows
+  """
+  Storey class
+  """
+  def __init__(self, name, surfaces, neighbours):
+    # todo: the information of the parent surface is lost -> need to recover it
+    self._name = name
+    self._surfaces = surfaces
+    self._thermal_boundaries = None
+    self._virtual_surfaces = None
+    self._thermal_zone = None
+    self._neighbours = neighbours
+
+  @property
+  def name(self):
+    """
+    Storey's name
+    :return: str
+    """
+    return self._name
+
+  @property
+  def surfaces(self) -> List[Surface]:
+    """
+    External surfaces enclosing the storey
+    :return: [Surface]
+    """
+    return self._surfaces
+
+  @property
+  def neighbours(self):
+    """
+    Neighbour storeys' names
+    :return: [str]
+    """
+    return self._neighbours
+
+  @property
+  def thermal_boundaries(self) -> List[ThermalBoundary]:
+    """
+    Thermal boundaries bounding the thermal zone
+    :return: [ThermalBoundary]
+    """
+    # todo: it cannot be, it creates a loop between thermal boundaries and thermal zones
+    if self._thermal_boundaries is None:
+      self._thermal_boundaries = []
+      for surface in self.surfaces:
+        if surface.type != cte.INTERIOR_WALL or surface.type != cte.INTERIOR_SLAB:
+          # external thermal boundary -> only one thermal zone
+          delimits = [self.thermal_zone]
+        else:
+          # internal thermal boundary -> two thermal zones
+          grad = np.rad2deg(surface.inclination)
+          if grad >= 170:
+            delimits = [self.thermal_zone, self._neighbours[0]]
+          else:
+            delimits = [self._neighbours[1], self.thermal_zone]
+        self._thermal_boundaries.append(ThermalBoundary(surface, delimits))
+    return self._thermal_boundaries
+
+  @property
+  def virtual_surfaces(self) -> List[Surface]:
+    """
+    Internal surfaces enclosing the thermal zone
+    :return: [Surface]
+    """
+    if self._virtual_surfaces is None:
+      self._virtual_surfaces = []
+      for thermal_boundary in self.thermal_boundaries:
+        self._virtual_surfaces.append(thermal_boundary.virtual_internal_surface)
+    return self._virtual_surfaces
+
+  @property
+  def thermal_zone(self) -> ThermalZone:
+    """
+    Thermal zone inside the storey
+    :return: ThermalZone
+    """
+    if self._thermal_zone is None:
+      self._thermal_zone = ThermalZone(self.virtual_surfaces)
+    return self._thermal_zone
diff --git a/city_model_structure/attributes/surface.py b/city_model_structure/attributes/surface.py
index ab591249..9a1a035d 100644
--- a/city_model_structure/attributes/surface.py
+++ b/city_model_structure/attributes/surface.py
@@ -9,24 +9,27 @@ from __future__ import annotations
 import numpy as np
 import uuid
 from city_model_structure.attributes.polygon import Polygon
+from city_model_structure.attributes.plane import Plane
+from city_model_structure.attributes.point import Point
 from city_model_structure.attributes.pv_system import PvSystem
+from city_model_structure.attributes.thermal_boundary import ThermalBoundary
+import helpers.constants as cte
 
 
 class Surface:
   """
   Surface class
   """
-  def __init__(self, solid_polygon, perimeter_polygon, holes_polygons=None, surface_type=None, swr=None):
+  def __init__(self, solid_polygon, perimeter_polygon, holes_polygons=None, name=None, surface_type=None, swr=None,
+               is_child=False):
     self._type = surface_type
     self._swr = swr
-    self._name = None
+    self._name = name
     self._id = None
     self._azimuth = None
     self._inclination = None
     self._area_above_ground = None
     self._area_below_ground = None
-    # todo @Guille: what is parent?????????????
-    self._parent = None
     self._lower_corner = None
     self._upper_corner = None
     self._shared_surfaces = []
@@ -35,6 +38,8 @@ class Surface:
     self._holes_polygons = holes_polygons
     self._solid_polygon = solid_polygon
     self._pv_system_installed = None
+    self._inverse = None
+    self._thermal_boundary = None
 
   @property
   def name(self):
@@ -244,3 +249,43 @@ class Surface:
     :param value: PvSystem
     """
     self._pv_system_installed = value
+
+  @property
+  def inverse(self) -> Surface:
+    """
+    Returns the same surface pointing backwards
+    :return: Surface
+    """
+    if self._inverse is None:
+      new_solid_polygon = Polygon(self.solid_polygon.inverse)
+      new_perimeter_polygon = Polygon(self.perimeter_polygon.inverse)
+      new_holes_polygons = []
+      for hole in self.holes_polygons:
+        new_holes_polygons.append(Polygon(hole.inverse))
+      self._inverse = Surface(new_solid_polygon, new_perimeter_polygon, new_holes_polygons, cte.VIRTUAL_INTERNAL)
+    return self._inverse
+
+  @property
+  def associated_thermal_boundary(self) -> ThermalBoundary:
+    """
+    Thermal boundary associated to this surface considered as the external face
+    :return: ThermalBoundary
+    """
+    if self._thermal_boundary is None:
+      self._thermal_boundary = ThermalBoundary(self, delimits)
+    return self._thermal_boundary
+
+  def shared_surfaces(self):
+    # todo: check https://trimsh.org/trimesh.collision.html as an option to implement this method
+    raise NotImplementedError
+
+  def divide(self, z):
+    # todo: recheck this method for LoD3 (windows)
+    origin = Point([0, 0, z])
+    normal = np.array([0, 0, 1])
+    plane = Plane(normal=normal, origin=origin)
+    polygon = self.perimeter_polygon
+    part_1, part_2, intersection = polygon.divide(plane)
+    surface_child = Surface(part_1, part_1, name=self.name, surface_type=self.type, is_child=True)
+    rest_surface = Surface(part_2, part_2, name=self.name, surface_type=self.type, is_child=True)
+    return surface_child, rest_surface, intersection
diff --git a/city_model_structure/attributes/thermal_boundary.py b/city_model_structure/attributes/thermal_boundary.py
index ae1d8dfc..b02b97c5 100644
--- a/city_model_structure/attributes/thermal_boundary.py
+++ b/city_model_structure/attributes/thermal_boundary.py
@@ -4,13 +4,13 @@ SPDX - License - Identifier: LGPL - 3.0 - or -later
 Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca
 Contributors Pilar Monsalvete Alvarez de Uribarri pilar.monsalvete@concordia.ca
 """
-from typing import List, TypeVar
+from typing import List, TypeVar, Union
 from city_model_structure.attributes.layer import Layer
 from city_model_structure.attributes.thermal_opening import ThermalOpening
-from city_model_structure.attributes.surface import Surface
 
 ThermalZone = TypeVar('ThermalZone')
 Polygon = TypeVar('Polygon')
+Surface = TypeVar('Surface')
 
 
 class ThermalBoundary:
@@ -34,6 +34,8 @@ class ThermalBoundary:
     self._window_ratio = None
     self._refurbishment_measure = None
     self._surface_geometry = None
+    self._thickness = None
+    self._virtual_internal_surface = None
 
   @property
   def surface(self) -> Surface:
@@ -76,22 +78,18 @@ class ThermalBoundary:
     """
     return self._surface.solid_polygon.area
 
-  # todo: shouldn't be these two attributes come from the associated surface???
   @property
-  def area_above_ground(self):
+  def thickness(self):
     """
-    Thermal boundary area above ground in square meters
+    Thermal boundary thickness in meters
     :return: float
     """
-    return self._surface.area_above_ground
-
-  @property
-  def area_below_ground(self):
-    """
-    Thermal boundary area below ground in square meters
-    :return: float
-    """
-    return self._surface.area_below_ground
+    if self._thickness is None:
+      self._thickness = 0.0
+      if self.layers is not None:
+        for layer in self.layers:
+          self._thickness += layer.thickness
+    return self._thickness
 
   @property
   def outside_solar_absorptance(self):
@@ -295,9 +293,19 @@ class ThermalBoundary:
     self._he = value
 
   @property
-  def surface_geometry(self) -> Polygon:
+  def surface_geometry(self) -> Union[NotImplementedError, Polygon]:
     """
     Get the polygon that defines the thermal boundary
     :return: Polygon
     """
-    return self._surface_geometry
+    raise NotImplementedError
+
+  @property
+  def virtual_internal_surface(self) -> Surface:
+    """
+    Get the internal surface of the thermal boundary
+    :return: Surface
+    """
+    if self._virtual_internal_surface is None:
+      self._virtual_internal_surface = self.surface.inverse
+    return self._virtual_internal_surface
diff --git a/city_model_structure/building.py b/city_model_structure/building.py
index 13dfccd6..68fc6ca8 100644
--- a/city_model_structure/building.py
+++ b/city_model_structure/building.py
@@ -11,12 +11,13 @@ from typing import List
 import numpy as np
 import math
 from city_model_structure.attributes.surface import Surface
-from city_model_structure.attributes.thermal_boundary import ThermalBoundary
 from city_model_structure.attributes.thermal_zone import ThermalZone
 from city_model_structure.attributes.usage_zone import UsageZone
+from city_model_structure.attributes.storey import Storey
+from city_model_structure.attributes.polygon import Polygon
+from city_model_structure.attributes.point import Point
 from city_model_structure.city_object import CityObject
-from helpers.geometry_helper import GeometryHelper as gh
-from trimesh import Trimesh
+from helpers import constants as cte
 
 
 class Building(CityObject):
@@ -24,7 +25,7 @@ class Building(CityObject):
   Building(CityObject) class
   """
   def __init__(self, name, lod, surfaces, year_of_construction, function,
-               city_lower_corner, terrains=None, zones_surfaces_ids=None):
+               city_lower_corner, terrains=None):
     super().__init__(name, lod, surfaces, city_lower_corner)
     self._basement_heated = None
     self._attic_heated = None
@@ -35,25 +36,12 @@ class Building(CityObject):
     self._storeys_above_ground = None
     self._floor_area = None
     self._roof_type = None
+    self._thermal_zones = []
     self._usage_zones = []
     self._type = 'building'
     self._heating = dict()
     self._cooling = dict()
     self._eave_height = None
-
-    self._thermal_zones = []
-    if zones_surfaces_ids is not None:
-      for zone_surfaces_ids in zones_surfaces_ids:
-        zone_surfaces = []
-        for surface_id in zone_surfaces_ids:
-          zone_surfaces.append(self.surface(surface_id))
-        self._thermal_zones.append(ThermalZone(zone_surfaces))
-    else:
-      zone_surfaces = surfaces
-      self._thermal_zones.append(ThermalZone(zone_surfaces))
-    for t_zones in self._thermal_zones:
-      t_zones.bounded = [ThermalBoundary(s, [t_zones]) for s in t_zones.surfaces]
-
     self._grounds = []
     self._roofs = []
     self._walls = []
@@ -311,40 +299,79 @@ class Building(CityObject):
     return self._eave_height
 
   @property
-  def storeys(self) -> [Trimesh]:
+  def storeys(self) -> [Storey]:
     """
-    subsections of building trimesh by storage in case of no interiors defined
-    :return: [Trimesh]
+    subsections of building trimesh by storey in case of no interiors defined
+    :return: [Storey]
     """
-    trimesh = self.simplified_polyhedron.trimesh
-    if self.average_storey_height is None:
-      if self.storeys_above_ground is None or self.storeys_above_ground <= 0:
-        sys.stderr.write('Warning: not enough information to divide building into storeys, '
-                         'either number of storeys or average storey height must be provided.\n')
-        return [trimesh]
-      else:
-        number_of_storeys = int(self.storeys_above_ground)
-      height = self.eave_height / number_of_storeys
-    else:
-      height = self.average_storey_height
-      if self.storeys_above_ground is not None:
-        number_of_storeys = int(self.storeys_above_ground)
-      else:
-        number_of_storeys = math.floor(float(self.eave_height) / height) + 1
-        last_storey_height = float(self.eave_height) - height*(number_of_storeys-1)
-        if last_storey_height < 0.3*height:
-          number_of_storeys -= 1
+    number_of_storeys, height = self._calculate_number_storeys_and_height(self.average_storey_height, self.eave_height,
+                                                                          self.storeys_above_ground)
+    if number_of_storeys == 0:
+      return Storey('storey_0', self.surfaces, [None, None])
 
     storeys = []
-    for n in range(0, number_of_storeys - 1):
-      point_plane = [self.city_object_lower_corner[0], self.city_object_lower_corner[1],
-                     self.city_object_lower_corner[2] + height * (n + 1)]
-      normal = [0, 0, -1]
-      storey, trimesh = gh.divide_mesh_by_plane(trimesh, normal, point_plane)
-      storeys.append(storey)
-    storeys.append(trimesh)
+    surfaces_child_last_storey = []
+    rest_surfaces = []
+    for i in range(0, number_of_storeys-1):
+      name = 'storey_' + str(i)
+      surfaces_child = []
+      if i == 0:
+        neighbours = [None, 'storey_1']
+        for surface in self.surfaces:
+          if surface.type == cte.GROUND:
+            surfaces_child.append(surface)
+          else:
+            rest_surfaces.append(surface)
+      else:
+        neighbours = ['storey_' + str(i-1), 'storey_' + str(i+1)]
+      height_division = self.lower_corner[2] + height*(i+1)
+      intersections = []
+      for surface in rest_surfaces:
+        if surface.type == cte.ROOF:
+          if height_division >= surface.upper_corner[2] > height_division-height:
+            surfaces_child.append(surface)
+          else:
+            surfaces_child_last_storey.append(surface)
+        else:
+          surface_child, rest_surface, intersection = surface.divide(height_division)
+          surfaces_child.append(surface_child)
+          intersections.extend(intersection)
+          if i == number_of_storeys-2:
+            surfaces_child_last_storey.append(rest_surface)
+      points = []
+      for intersection in intersections:
+        points.append(intersection[1])
+      coordinates = self._intersections_to_coordinates(intersections)
+      polygon = Polygon(coordinates)
+      ceiling = Surface(polygon, polygon, surface_type=cte.INTERIOR_SLAB)
+      surfaces_child.append(ceiling)
+      storeys.append(Storey(name, surfaces_child, neighbours))
+    name = 'storey_' + str(number_of_storeys-1)
+    neighbours = ['storey_' + str(number_of_storeys-2), None]
+    storeys.append(Storey(name, surfaces_child_last_storey, neighbours))
     return storeys
 
+  @staticmethod
+  def _calculate_number_storeys_and_height(average_storey_height, eave_height, storeys_above_ground):
+    if average_storey_height is None:
+      if storeys_above_ground is None or storeys_above_ground <= 0:
+        sys.stderr.write('Warning: not enough information to divide building into storeys, '
+                         'either number of storeys or average storey height must be provided.\n')
+        return 0, 0
+      else:
+        number_of_storeys = int(storeys_above_ground)
+      height = eave_height / number_of_storeys
+    else:
+      height = average_storey_height
+      if storeys_above_ground is not None:
+        number_of_storeys = int(storeys_above_ground)
+      else:
+        number_of_storeys = math.floor(float(eave_height) / height) + 1
+        last_storey_height = float(eave_height) - height*(number_of_storeys-1)
+        if last_storey_height < 0.3*height:
+          number_of_storeys -= 1
+    return number_of_storeys, height
+
   @property
   def roof_type(self):
     """
@@ -379,3 +406,32 @@ class Building(CityObject):
   @pv_plus_hp_installation.setter
   def pv_plus_hp_installation(self, value):
     self._pv_plus_hp_installation = value
+
+  @staticmethod
+  def _intersections_to_coordinates(edges_list):
+    # todo: this method is horrible, the while loop needs to be improved
+    points = [Point(edges_list[0][0]), Point(edges_list[0][1])]
+    found_edges = []
+    j = 0
+    while j < len(points)-1:
+      for i in range(1, len(edges_list)):
+        if i not in found_edges:
+          point_2 = points[len(points) - 1]
+          point_1 = Point(edges_list[i][0])
+          found = False
+          if point_1.distance_to_point(point_2) <= 1e-10:
+            points.append(Point(edges_list[i][1]))
+            found_edges.append(i)
+            found = True
+          if not found:
+            point_1 = Point(edges_list[i][1])
+            if point_1.distance_to_point(point_2) <= 1e-10:
+              points.append(Point(edges_list[i][0]))
+              found_edges.append(i)
+      j += 1
+
+    points.remove(points[len(points)-1])
+    array_points = []
+    for point in points:
+      array_points.append(point.coordinates)
+    return np.array(array_points)
diff --git a/city_model_structure/city_object.py b/city_model_structure/city_object.py
index 6009a5c8..aa08dedc 100644
--- a/city_model_structure/city_object.py
+++ b/city_model_structure/city_object.py
@@ -202,7 +202,7 @@ class CityObject:
     self._beam = value
 
   @property
-  def city_object_lower_corner(self):
+  def lower_corner(self):
     """
     City object lower corner coordinates [x, y, z]
     """
diff --git a/helpers/constants.py b/helpers/constants.py
index 31597435..bfbb3bdf 100644
--- a/helpers/constants.py
+++ b/helpers/constants.py
@@ -18,6 +18,7 @@ ATTIC_FLOOR = 'Attic floor'
 ROOF = 'Roof'
 INTERIOR_SLAB = 'Interior slab'
 INTERIOR_WALL = 'Interior wall'
+VIRTUAL_INTERNAL = 'Virtual internal'
 WINDOW = 'Window'
 DOOR = 'Door'
 SKYLIGHT = 'Skylight'
diff --git a/non_functional_tests/test_construction_factory.py b/non_functional_tests/test_construction_factory.py
index e0613e56..8b8a3432 100644
--- a/non_functional_tests/test_construction_factory.py
+++ b/non_functional_tests/test_construction_factory.py
@@ -51,6 +51,7 @@ class TestConstructionFactory(TestCase):
     for building in city.buildings:
       self.assertIsNotNone(building.average_storey_height, 'average_storey_height is none')
       self.assertIsNotNone(building.storeys_above_ground, 'storeys_above_ground is none')
+      self.assertIsNot(len(building.thermal_zones), 0, 'no building thermal_zones defined')
       for thermal_zone in building.thermal_zones:
         self.assertIsNotNone(thermal_zone.effective_thermal_capacity, 'effective_thermal_capacity is none')
         self.assertIsNotNone(thermal_zone.additional_thermal_bridge_u_value,
@@ -81,6 +82,7 @@ class TestConstructionFactory(TestCase):
     for building in city.buildings:
       self.assertIsNotNone(building.average_storey_height, 'average_storey_height is none')
       self.assertIsNotNone(building.storeys_above_ground, 'storeys_above_ground is none')
+      self.assertIsNot(len(building.thermal_zones), 0, 'no building thermal_zones defined')
       for thermal_zone in building.thermal_zones:
         self.assertIsNotNone(thermal_zone.effective_thermal_capacity, 'effective_thermal_capacity is none')
         self.assertIsNotNone(thermal_zone.additional_thermal_bridge_u_value,
diff --git a/non_functional_tests/test_geometry_factory.py b/non_functional_tests/test_geometry_factory.py
index 62e30a64..5f283d84 100644
--- a/non_functional_tests/test_geometry_factory.py
+++ b/non_functional_tests/test_geometry_factory.py
@@ -71,7 +71,7 @@ class TestGeometryFactory(TestCase):
       self.assertIsNotNone(building.usage_zones, 'building usage_zones is none')
       self.assertIsNone(building.average_storey_height, 'building average_storey_height is not none')
       self.assertIsNone(building.storeys_above_ground, 'building storeys_above_ground is not none')
-      self.assertIsNotNone(building.thermal_zones, 'building thermal_zones is none')
+      self.assertIsNot(len(building.thermal_zones), 0, 'no building thermal_zones defined')
       self.assertIsNotNone(building.type, 'building type is none')
       self.assertIsNotNone(building.max_height, 'building max_height is none')
       self.assertIsNotNone(building.floor_area, 'building floor_area is none')
@@ -108,6 +108,7 @@ class TestGeometryFactory(TestCase):
     file = 'one_building_in_kelowna.gml'
     city = self._get_citygml(file)
     for building in city.buildings:
+      self.assertIsNot(len(building.thermal_zones), 0, 'no building thermal_zones defined')
       for thermal_zone in building.thermal_zones:
         self.assertIsNotNone(thermal_zone.surfaces, 'thermal_zone surfaces is none')
         self.assertIsNotNone(thermal_zone.bounded, 'thermal_zone bounded is none')
@@ -135,7 +136,9 @@ class TestGeometryFactory(TestCase):
     file = 'one_building_in_kelowna.gml'
     city = self._get_citygml(file)
     for building in city.buildings:
+      self.assertIsNot(len(building.thermal_zones), 0, 'no building thermal_zones defined')
       for thermal_zone in building.thermal_zones:
+        self.assertIsNot(len(thermal_zone.bounded), 0, 'no building thermal_boundaries defined')
         for thermal_boundary in thermal_zone.bounded:
           self.assertIsNotNone(thermal_boundary.surface, 'thermal_boundary surface is none')
           self.assertIsNotNone(thermal_boundary.type, 'thermal_boundary type is none')
@@ -154,7 +157,9 @@ class TestGeometryFactory(TestCase):
     file = 'one_building_in_kelowna.gml'
     city = self._get_citygml(file)
     for building in city.buildings:
+      self.assertIsNot(len(building.thermal_zones), 0, 'no building thermal_zones defined')
       for thermal_zone in building.thermal_zones:
+        self.assertIsNot(len(thermal_zone.bounded), 0, 'no building thermal_boundaries defined')
         for thermal_boundary in thermal_zone.bounded:
           for thermal_opening in thermal_boundary.thermal_openings:
             self.assertIsNone(thermal_opening.frame_ratio, 'thermal_opening frame_ratio was initialized')
@@ -200,6 +205,31 @@ class TestGeometryFactory(TestCase):
     for building in city.buildings:
       self.assertRaises(Exception, lambda: self._internal_function(function_format, building.function))
 
+  def test_citygml_storeys(self):
+    """
+    Test division by storeys of buildings
+    :return: None
+    """
+
+    file = 'one_building_in_kelowna.gml'
+    city = self._get_citygml(file)
+    for building in city.buildings:
+      print('building')
+      for surface in building.surfaces:
+        print(surface.name)
+        print(surface.type)
+        print(surface.perimeter_polygon.area)
+      building.average_storey_height = 1.5
+      building.storeys_above_ground = 2
+      storeys = building.storeys
+      for storey in storeys:
+        print(storey.name)
+        print(storey.neighbours)
+        for surface in storey.surfaces:
+          print(surface.name)
+          print(surface.type)
+          print(surface.perimeter_polygon.area)
+
   # obj
   def test_import_obj(self):
     file = 'kelowna.obj'
diff --git a/non_functional_tests/test_schedules_factory.py b/non_functional_tests/test_schedules_factory.py
index d50d7f1e..f04dbd89 100644
--- a/non_functional_tests/test_schedules_factory.py
+++ b/non_functional_tests/test_schedules_factory.py
@@ -40,5 +40,6 @@ class TestSchedulesFactory(TestCase):
     occupancy_handler = 'comnet'
     SchedulesFactory(occupancy_handler, city).enrich()
     for building in city.buildings:
+      self.assertIsNot(len(building.usage_zones), 0, 'no building usage_zones defined')
       for usage_zone in building.usage_zones:
         self.assertTrue(usage_zone.schedules)
diff --git a/non_functional_tests/test_usage_factory.py b/non_functional_tests/test_usage_factory.py
index 2d13d4e5..d48b04d6 100644
--- a/non_functional_tests/test_usage_factory.py
+++ b/non_functional_tests/test_usage_factory.py
@@ -42,7 +42,7 @@ class TestUsageFactory(TestCase):
     # case 1: HFT
     UsageFactory('hft', city).enrich()
     for building in city.buildings:
-      self.assertIsNotNone(building.usage_zones, 'usage_zones not created')
+      self.assertIsNot(len(building.usage_zones), 0, 'no building usage_zones defined')
       for usage_zone in building.usage_zones:
         self.assertIsNotNone(usage_zone.usage, 'usage is none')
         self.assertIsNotNone(usage_zone.internal_gains, 'usage is none')