diff --git a/city_model_structure/city.py b/city_model_structure/city.py
new file mode 100644
index 00000000..c7a169c5
--- /dev/null
+++ b/city_model_structure/city.py
@@ -0,0 +1,43 @@
+from city_model_structure.city_object import CityObject
+from typing import List, Union
+
+
+class City:
+  def __init__(self, lower_corner, upper_corner, city_objects, srs_name):
+    self._city_objects = None
+    self._name = None
+    self._lower_corner = lower_corner
+    self._upper_corner = upper_corner
+    self._city_objects = city_objects
+    self._srs_name = srs_name
+
+  @property
+  def city_objects(self) -> List[CityObject]:
+    return self._city_objects
+
+  @property
+  def lower_corner(self):
+    return self._lower_corner
+
+  @property
+  def upper_corner(self):
+    return self._upper_corner
+
+  def city_object(self, name) -> Union[CityObject, None]:
+    for c in self.city_objects:
+      if c.name == name:
+        return c
+    return None
+
+  @property
+  def srs_name(self):
+    return self._srs_name
+
+  @property
+  def name(self):
+    return self._name
+
+  @name.setter
+  def name(self, value):
+    self._name = value
+
diff --git a/city_model_structure/city_object.py b/city_model_structure/city_object.py
new file mode 100644
index 00000000..67af8767
--- /dev/null
+++ b/city_model_structure/city_object.py
@@ -0,0 +1,193 @@
+from matplotlib import pylab
+from city_model_structure.polyhedron import Polyhedron
+from city_model_structure.thermal_zone import ThermalZone
+from city_model_structure.thermal_boundary import ThermalBoundary
+from city_model_structure.surface import Surface
+from shapely import ops
+from shapely.geometry import MultiPolygon
+import numpy as np
+from pathlib import Path
+import matplotlib.patches as patches
+from helpers.geometry import Geometry
+from city_model_structure.usage_zone import UsageZone
+from typing import Union, List
+
+
+class CityObject:
+  def __init__(self, name, lod, surfaces, terrains, year_of_construction, function, attic_heated=0, basement_heated=0):
+    self._name = name
+    self._lod = lod
+    self._surfaces = surfaces
+    self._polyhedron = None
+    self._basement_heated = basement_heated
+    self._attic_heated = attic_heated
+    self._terrains = terrains
+    self._year_of_construction = year_of_construction
+    self._function = function
+    self._geometry = Geometry()
+    self._average_storey_height = None
+    self._storeys_above_ground = None
+    self._foot_print = None
+    self._usage_zones = []
+
+    # ToDo: Check this for LOD4
+    self._thermal_zones = []
+    if self.lod < 8:
+      # for lod under 4 is just one thermal zone
+      self._thermal_zones.append(ThermalZone(self.surfaces))
+    for t in self._thermal_zones:
+      t.bounded = [ThermalBoundary(s, [t]) for s in t.surfaces]
+    surface_id = 0
+    for s in self._surfaces:
+      s.parent(self, surface_id)
+      surface_id += 1
+
+  @property
+  def usage_zones(self) -> List[UsageZone]:
+    return self._usage_zones
+
+  @usage_zones.setter
+  def usage_zones(self, value):
+    self._usage_zones = value
+
+  @property
+  def terrains(self) -> List[Surface]:
+    return self._terrains
+
+  @property
+  def attic_heated(self):
+    return self._attic_heated
+
+  @attic_heated.setter
+  def attic_heated(self, value):
+    self._attic_heated = value
+
+  @property
+  def basement_heated(self):
+    return self._basement_heated
+
+  @basement_heated.setter
+  def basement_heated(self, value):
+    self._attic_heated = value
+
+  @property
+  def name(self):
+    return self._name
+
+  @property
+  def lod(self):
+    return self._lod
+
+  @property
+  def surfaces(self) -> List[Surface]:
+    return self._surfaces
+
+  def surface(self, name) -> Union[Surface, None]:
+    for s in self.surfaces:
+      if s.name == name:
+        return s
+    return None
+
+  @property
+  def thermal_zones(self) -> List[ThermalZone]:
+    return self._thermal_zones
+
+  @property
+  def volume(self):
+    if self._polyhedron is None:
+      self._polyhedron = Polyhedron(self.surfaces)
+    return self._polyhedron.volume
+
+  @property
+  def heated_volume(self):
+    if self._polyhedron is None:
+      self._polyhedron = Polyhedron(self.surfaces)
+    # ToDo: this need to be the calculated based on the basement and attic heated values
+    return self._polyhedron.volume
+
+  def stl_export(self, path):
+    if self._polyhedron is None:
+      self._polyhedron = Polyhedron(self.surfaces)
+    full_path = (Path(path) / (self._name + '.stl')).resolve()
+    self._polyhedron.save(full_path)
+
+  @property
+  def year_of_construction(self):
+    return self._year_of_construction
+
+  @property
+  def function(self):
+    return self._function
+
+  @property
+  def average_storey_height(self):
+    return self._average_storey_height
+
+  @average_storey_height.setter
+  def average_storey_height(self, value):
+    self._average_storey_height = value
+
+  @property
+  def storeys_above_ground(self):
+    return self._storeys_above_ground
+
+  @storeys_above_ground.setter
+  def storeys_above_ground(self, value):
+    self._storeys_above_ground = value
+
+  @staticmethod
+  def _tuple_to_point(xy_tuple):
+    return [xy_tuple[0], xy_tuple[1], 0.0]
+
+  def _plot(self, polygon):
+    points = ()
+    for point_tuple in polygon.exterior.coords:
+      almost_equal = False
+      for point in points:
+        p1 = CityObject._tuple_to_point(point)
+        p2 = CityObject._tuple_to_point(point_tuple)
+        if self._geometry.almost_equal(p1, p2):
+          almost_equal = True
+          break
+      if not almost_equal:
+        points = points + (point_tuple,)
+    points = points + (points[0],)
+    pylab.scatter([point[0] for point in points], [point[1] for point in points])
+    pylab.gca().add_patch(patches.Polygon(points, closed=True, fill=True))
+    pylab.grid()
+    pylab.show()
+
+  @property
+  def foot_print(self) -> Surface:
+    if self._foot_print is None:
+      shapelys = []
+      union = None
+      for surface in self.surfaces:
+        if surface.shapely.is_empty or not surface.shapely.is_valid:
+          continue
+        shapelys.append(surface.shapely)
+        union = ops.unary_union(shapelys)
+        shapelys = [union]
+      if type(union) == MultiPolygon:
+        Exception('foot print returns a multipolygon')
+      points_list = []
+      for point_tuple in union.exterior.coords:
+        # ToDo: should be Z 0.0 or min Z?
+        point = CityObject._tuple_to_point(point_tuple)
+        almost_equal = False
+        for existing_point in points_list:
+          if self._geometry.almost_equal(point, existing_point):
+            almost_equal = True
+            break
+        if not almost_equal:
+          points_list.append(point)
+      points_list = np.reshape(points_list, len(points_list) * 3)
+      points = np.array_str(points_list).replace('[', '').replace(']', '')
+      self._foot_print = Surface(points, remove_last=False, is_projected=True)
+    return self._foot_print
+
+  @property
+  def max_height(self):
+    if self._polyhedron is None:
+      self._polyhedron = Polyhedron(self.surfaces)
+    return self._polyhedron.max_z
diff --git a/city_model_structure/internal_gains.py b/city_model_structure/internal_gains.py
new file mode 100644
index 00000000..3f86d7c8
--- /dev/null
+++ b/city_model_structure/internal_gains.py
@@ -0,0 +1,38 @@
+class InternalGains:
+  def __init__(self):
+    self._average_internal_gain_w_m2 = None
+    self._convective_fraction = None
+    self._radiative_fraction = None
+    self._latent_fraction = None
+
+  @property
+  def average_internal_gain_w_m2(self):
+    return self._average_internal_gain_w_m2
+
+  @average_internal_gain_w_m2.setter
+  def average_internal_gain_w_m2(self, value):
+    self._average_internal_gain_w_m2 = value
+
+  @property
+  def convective_fraction(self):
+    return self._convective_fraction
+
+  @convective_fraction.setter
+  def convective_fraction(self, value):
+    self._convective_fraction = value
+
+  @property
+  def radiative_fraction(self):
+    return self._radiative_fraction
+
+  @radiative_fraction.setter
+  def radiative_fraction(self, value):
+    self._radiative_fraction = value
+
+  @property
+  def latent_fraction(self):
+    return self._latent_fraction
+
+  @latent_fraction.setter
+  def latent_fraction(self, value):
+    self._latent_fraction = value
diff --git a/city_model_structure/layer.py b/city_model_structure/layer.py
new file mode 100644
index 00000000..bb96b161
--- /dev/null
+++ b/city_model_structure/layer.py
@@ -0,0 +1,23 @@
+from city_model_structure.material import Material
+
+
+class Layer:
+  def __init__(self):
+    self._material = None
+    self._thickness_m = None
+
+  @property
+  def material(self) -> Material:
+    return self._material
+
+  @material.setter
+  def material(self, value):
+    self._material = value
+
+  @property
+  def thickness_m(self):
+    return self._thickness_m
+
+  @thickness_m.setter
+  def thickness_m(self, value):
+    self._thickness_m = value
diff --git a/city_model_structure/material.py b/city_model_structure/material.py
new file mode 100644
index 00000000..2473c9ad
--- /dev/null
+++ b/city_model_structure/material.py
@@ -0,0 +1,75 @@
+class Material:
+  def __init__(self):
+    # ToDo: construct this class
+    self._conductivity_wm_k = None
+    self._specific_heat_jkg_k = None
+    self._density_kg_m3 = None
+    self._solar_absorptance = None
+    self._thermal_absorptance = None
+    self._visible_absorptance = None
+    self._no_mass = False
+    self._thermal_resistance_m2k_w = None
+
+  @property
+  def conductivity_wm_k(self):
+    return self._conductivity_wm_k
+
+  @conductivity_wm_k.setter
+  def conductivity_wm_k(self, value):
+    self._conductivity_wm_k = value
+
+  @property
+  def specific_heat_jkg_k(self):
+    return self._specific_heat_jkg_k
+
+  @specific_heat_jkg_k.setter
+  def specific_heat_jkg_k(self, value):
+    self._specific_heat_jkg_k = value
+
+  @property
+  def density_kg_m3(self):
+    return self._density_kg_m3
+
+  @density_kg_m3.setter
+  def density_kg_m3(self, value):
+    self._density_kg_m3 = value
+
+  @property
+  def solar_absorptance(self):
+    return self._solar_absorptance
+
+  @solar_absorptance.setter
+  def solar_absorptance(self, value):
+    self._solar_absorptance = value
+
+  @property
+  def thermal_absorptance(self):
+    return self._thermal_absorptance
+
+  @thermal_absorptance.setter
+  def thermal_absorptance(self, value):
+    self._thermal_absorptance = value
+
+  @property
+  def visible_absorptance(self):
+    return self._visible_absorptance
+
+  @visible_absorptance.setter
+  def visible_absorptance(self, value):
+    self._visible_absorptance = value
+
+  @property
+  def no_mass(self):
+    return self._no_mass
+
+  @no_mass.setter
+  def no_mass(self, value):
+    self._no_mass = value
+
+  @property
+  def thermal_resistance_m2k_w(self):
+    return self._thermal_resistance_m2k_w
+
+  @thermal_resistance_m2k_w.setter
+  def thermal_resistance_m2k_w(self, value):
+    self._thermal_resistance_m2k_w = value
diff --git a/city_model_structure/polyhedron.py b/city_model_structure/polyhedron.py
new file mode 100644
index 00000000..24385791
--- /dev/null
+++ b/city_model_structure/polyhedron.py
@@ -0,0 +1,74 @@
+import numpy as np
+import stl
+from helpers.geometry import Geometry
+
+
+class Polyhedron:
+  def __init__(self, surfaces):
+    self._surfaces = [s for s in surfaces]
+    self._polygons = [s.polygon for s in surfaces]
+    self._polyhedron = None
+    self._volume = None
+    self._faces = None
+    self._vertices = None
+    self._mesh = None
+    self._geometry = Geometry()
+
+  def _position_of(self, point):
+    vertices = self.vertices
+    for i in range(len(vertices)):
+      if self._geometry.almost_equal(vertices[i], point):
+        return i
+    return -1
+
+  @property
+  def vertices(self):
+    if self._vertices is None:
+      vertices, self._vertices = [], []
+      [vertices.extend(s.points) for s in self._surfaces]
+      for v1 in vertices:
+        found = False
+        for v2 in self._vertices:
+          found = False
+          if self._geometry.almost_equal(v1, v2):
+            found = True
+            break
+        if not found:
+          self._vertices.append(v1)
+      self._vertices = np.asarray(self._vertices)
+    return self._vertices
+
+  @property
+  def faces(self):
+    if self._faces is None:
+      self._faces = []
+      for s in self._surfaces:
+        face = []
+        points = s.points
+        for p in points:
+          face.append(self._position_of(p))
+        self._faces.append(face)
+      self._faces = np.asarray(self._faces)
+    return self._faces
+
+  @property
+  def _polyhedron_mesh(self):
+    if self._mesh is None:
+      self._mesh = stl.mesh.Mesh(np.zeros(self.faces.shape[0], dtype=stl.mesh.Mesh.dtype))
+      for i, f in enumerate(self.faces):
+        for j in range(3):
+          self._mesh.vectors[i][j] = self.vertices[f[j], :]
+    return self._mesh
+
+  @property
+  def volume(self):
+    if self._volume is None:
+      self._volume, cog, inertia = self._polyhedron_mesh.get_mass_properties()
+    return self._volume
+
+  @property
+  def max_z(self):
+    return self._polyhedron_mesh.z.max()
+
+  def save(self, full_path):
+    self._polyhedron_mesh.save(full_path)
diff --git a/city_model_structure/surface.py b/city_model_structure/surface.py
new file mode 100644
index 00000000..d3e8cbf8
--- /dev/null
+++ b/city_model_structure/surface.py
@@ -0,0 +1,176 @@
+from __future__ import annotations
+import numpy as np
+import pyny3d.geoms as pn
+from helpers.geometry import Geometry
+
+
+class Surface:
+  def __init__(self, coordinates, surface_type=None, name=None, swr='0.2', remove_last=True, is_projected=False):
+    self._coordinates = coordinates
+    self._type = surface_type
+    self._name = name
+    self._swr = swr
+    self._remove_last = remove_last
+    self._is_projected = is_projected
+    self._geometry = Geometry()
+    self._polygon = None
+    self._area = None
+    self._points = None
+    self._points_list = None
+    self._normal = None
+    self._azimuth = None
+    self._inclination = None
+    self._area_above_ground = None
+    self._area_below_ground = None
+    self._parent = None
+    self._shapely = None
+    self._projected_surface = None
+    self._global_irradiance_hour = np.zeros(8760)
+    self._global_irradiance_month = np.zeros(12)
+
+  def parent(self, parent, surface_id):
+    self._parent = parent
+    self._name = str(surface_id)
+
+  @property
+  def name(self):
+    if self._name is None:
+      raise Exception('surface has no name')
+    return self._name
+
+  @property
+  def swr(self):
+    return self._swr
+
+  @swr.setter
+  def swr(self, value):
+    self._swr = value
+
+  @property
+  def points(self):
+    if self._points is None:
+      self._points = np.fromstring(self._coordinates, dtype=float, sep=' ')
+      self._points = Geometry.to_points_matrix(self._points, self._remove_last)
+    return self._points
+
+  @property
+  def points_list(self):
+    if self._points_list is None:
+      s = self.points
+      self._points_list = np.reshape(s, len(s) * 3)
+    return self._points_list
+
+  @property
+  def polygon(self):
+    if self._polygon is None:
+      try:
+        self._polygon = pn.Polygon(self.points)
+      except ValueError:
+        # is not really a polygon but a line so just return none
+        self._polygon = None
+    return self._polygon
+
+  @property
+  def area(self):
+    if self._area is None:
+      self._area = self.polygon.get_area()
+    return self._area
+
+  def is_almost_same_terrain(self, terrain_points, ground_points):
+    equal = 0
+    for t in terrain_points:
+      for g in ground_points:
+        if self._geometry.almost_equal(t, g):
+          equal += 1
+    return equal == len(terrain_points)
+
+  @property
+  def is_terrain(self):
+    for t_points in self._parent.terrains:
+      if len(t_points) == len(self.points) and self.is_almost_same_terrain(t_points, self.points):
+        return True
+    return False
+
+  @property
+  def area_above_ground(self):
+    if self._area_above_ground is None:
+      self._area_above_ground = self.area - self.area_below_ground
+    return self._area_above_ground
+
+  @property
+  def area_below_ground(self):
+    if self._area_below_ground is None:
+      self._area_below_ground = 0.0
+      if self.is_terrain:
+        self._area_below_ground = self.area
+    return self._area_below_ground
+
+  @property
+  def normal(self):
+    if self._normal is None:
+      points = self.points
+      n = np.cross(points[1] - points[0], points[2] - points[0])
+      self._normal = n / np.linalg.norm(n)
+    return self._normal
+
+  @property
+  def azimuth(self):
+    if self._azimuth is None:
+      normal = self.normal
+      self._azimuth = np.arctan2(normal[1], normal[0])
+    return self._azimuth
+
+  @property
+  def inclination(self):
+    if self._inclination is None:
+      self._inclination = np.arccos(self.normal[2])
+    return self._inclination
+
+  @property
+  def type(self):
+    if self._type is None:
+      grad = np.rad2deg(self.inclination)
+      if 170 <= grad:
+        self._type = 'Ground'
+      elif 80 <= grad <= 100:
+        self._type = 'Wall'
+      else:
+        self._type = 'Roof'
+    return self._type
+
+  @property
+  def global_irradiance_hour(self):
+    return self._global_irradiance_hour
+
+  @global_irradiance_hour.setter
+  def global_irradiance_hour(self, value):
+    self._global_irradiance_hour = value
+
+  @property
+  def global_irradiance_month(self):
+    return self._global_irradiance_month
+
+  @global_irradiance_month.setter
+  def global_irradiance_month(self, value):
+    self._global_irradiance_month = value
+
+  @property
+  def shapely(self):
+    if self.polygon is None:
+      return None
+    if self._shapely is None:
+      self._shapely = self.polygon.get_shapely()
+    return self._shapely
+
+  @property
+  def projection(self) -> Surface:
+    if self._is_projected:
+      return self
+    if self._projected_surface is None:
+      coordinates = ''
+      for coordinate in self.shapely.exterior.coords:
+        if coordinates != '':
+          coordinates = coordinates + ' '
+        coordinates = coordinates + str(coordinate[0]) + ' ' + str(coordinate[1]) + ' 0.0'
+      self._projected_surface = Surface(coordinates)
+    return self._projected_surface
diff --git a/city_model_structure/thermal_boundary.py b/city_model_structure/thermal_boundary.py
new file mode 100644
index 00000000..2c7e62fc
--- /dev/null
+++ b/city_model_structure/thermal_boundary.py
@@ -0,0 +1,132 @@
+from city_model_structure.thermal_opening import ThermalOpening
+from city_model_structure.surface import Surface
+from city_model_structure.layer import Layer
+import helpers.assumptions as assumptions
+from typing import List
+
+
+class ThermalBoundary:
+  def __init__(self, surface, delimits):
+    self._surface = surface
+    self._delimits = delimits
+    # ToDo: up to at least LOD2 will be just one thermal opening per Thermal boundary, review for LOD3 and LOD4
+    self._thermal_openings = [ThermalOpening()]
+    self._layers = None
+    self._outside_solar_absorptance = None
+    self._outside_thermal_absorptance = None
+    self._outside_visible_absorptance = None
+    self._window_ratio = None
+    self._u_value = None
+    self._window_area = None
+    self._shortwave_reflectance = None
+
+  @property
+  def delimits(self) -> List[Surface]:
+    return self._delimits
+
+  @property
+  def azimuth(self):
+    return self._surface.azimuth
+  
+  @property
+  def inclination(self):
+    return self._surface.inclination
+  
+  @property
+  def area(self):
+    return self._surface.area
+
+  @property
+  def area_above_ground(self):
+    return self._surface.area_above_ground
+
+  @property
+  def area_below_ground(self):
+    return self._surface.area_below_ground
+
+  @property
+  def outside_solar_absorptance(self):
+    return self._outside_solar_absorptance
+
+  @outside_solar_absorptance.setter
+  def outside_solar_absorptance(self, value):
+    self._outside_solar_absorptance = value
+
+  @property
+  def outside_thermal_absorptance(self):
+    return self._outside_thermal_absorptance
+
+  @outside_thermal_absorptance.setter
+  def outside_thermal_absorptance(self, value):
+    self._outside_thermal_absorptance = value
+
+  @property
+  def outside_visible_absorptance(self):
+    return self._outside_visible_absorptance
+
+  @outside_visible_absorptance.setter
+  def outside_visible_absorptance(self, value):
+    self._outside_visible_absorptance = value
+
+  @property
+  def thermal_openings(self) -> List[ThermalOpening]:
+    return self._thermal_openings
+
+  @thermal_openings.setter
+  def thermal_openings(self, value):
+    self._thermal_openings = value
+
+  @property
+  def layers(self) -> List[Layer]:
+    return self._layers
+
+  @layers.setter
+  def layers(self, value):
+    self._layers = value
+
+  @property
+  def type(self):
+    return self._surface.type
+
+  @property
+  def window_ratio(self):
+    return self._window_ratio
+
+  @window_ratio.setter
+  def window_ratio(self, value):
+    self._window_ratio = value
+
+  @property
+  def window_area(self):
+    if self._window_area is None:
+      try:
+        self._window_area = float(self._surface.area) * float(self.window_ratio)
+      except TypeError:
+        raise Exception('Window ratio is not defined or invalid surface area')
+    return self._window_area
+
+  @property
+  def u_value(self):
+    if self._u_value is None:
+      h_i = assumptions.h_i
+      h_e = assumptions.h_e
+      r_value = 1.0/h_i + 1.0/h_e
+      try:
+        for layer in self.layers:
+          if layer.material.no_mass:
+            r_value += float(layer.material.thermal_resistance_m2k_w)
+          else:
+            r_value = r_value + float(layer.material.conductivity_wm_k)/float(layer.thickness_m)
+        self._u_value = 1.0/r_value
+      except TypeError:
+        raise Exception('Constructions layers are not initialized')
+    return self._u_value
+
+  @property
+  def shortwave_reflectance(self):
+    if self._shortwave_reflectance is None:
+      try:
+        self._shortwave_reflectance = 1.0 - float(self.outside_solar_absorptance)
+      except TypeError:
+        raise Exception('Outside solar absorptance is not initialized')
+    return self._shortwave_reflectance
diff --git a/city_model_structure/thermal_opening.py b/city_model_structure/thermal_opening.py
new file mode 100644
index 00000000..f4072afc
--- /dev/null
+++ b/city_model_structure/thermal_opening.py
@@ -0,0 +1,90 @@
+import helpers.assumptions as assumptions
+
+
+class ThermalOpening:
+  def __init__(self):
+    self._area = None
+    self._openable_ratio = None
+    self._conductivity_w_mk = None
+    self._frame_ratio = assumptions.frame_ratio
+    self._g_value = None
+    self._thickness_m = None
+    self._inside_reflectance = None
+    self._outside_reflectance = None
+    self._u_value = None
+
+  @property
+  def area(self):
+    return self._area
+
+  @area.setter
+  def area(self, value):
+    self._area = value
+
+  @property
+  def openable_ratio(self):
+    raise Exception('Not implemented')
+
+  @openable_ratio.setter
+  def openable_ratio(self, value):
+    raise Exception('Not implemented')
+
+  @property
+  def conductivity_w_mk(self):
+    return self._conductivity_w_mk
+
+  @conductivity_w_mk.setter
+  def conductivity_w_mk(self, value):
+    self._conductivity_w_mk = value
+
+  @property
+  def frame_ratio(self):
+    return self._frame_ratio
+
+  @frame_ratio.setter
+  def frame_ratio(self, value):
+    self._frame_ratio = value
+
+  @property
+  def g_value(self):
+    return self._g_value
+
+  @g_value.setter
+  def g_value(self, value):
+    self._g_value = value
+
+  @property
+  def thickness_m(self):
+    return self._thickness_m
+
+  @thickness_m.setter
+  def thickness_m(self, value):
+    self._thickness_m = value
+
+  @property
+  def inside_reflectance(self):
+    return self._inside_reflectance
+
+  @inside_reflectance.setter
+  def inside_reflectance(self, value):
+    self._inside_reflectance = value
+
+  @property
+  def outside_reflectance(self):
+    return self._outside_reflectance
+
+  @outside_reflectance.setter
+  def outside_reflectance(self, value):
+    self._outside_reflectance = value
+
+  @property
+  def u_value(self):
+    if self._u_value is None:
+      try:
+        h_i = assumptions.h_i
+        h_e = assumptions.h_e
+        r_value = 1 / h_i + 1 / h_e + float(self.conductivity_w_mk) / float(self.thickness_m)
+        self._u_value = 1 / r_value
+      except TypeError:
+        Exception('Thermal opening is not initialized')
+    return self._u_value
diff --git a/city_model_structure/thermal_zone.py b/city_model_structure/thermal_zone.py
new file mode 100644
index 00000000..4f9533db
--- /dev/null
+++ b/city_model_structure/thermal_zone.py
@@ -0,0 +1,86 @@
+from typing import List
+from city_model_structure.thermal_boundary import ThermalBoundary
+
+
+class ThermalZone:
+  def __init__(self, surfaces, heated=True, cooled=True):
+    self._surfaces = surfaces
+    self._floor_area = None
+    self._bounded = None
+    self._heated = heated
+    self._cooled = cooled
+    self._additional_thermal_bridge_u_value = None
+    self._effective_thermal_capacity = None
+    self._indirectly_heated_area_ratio = None
+    self._infiltration_rate_system_on = None
+    self._infiltration_rate_system_off = None
+
+  @property
+  def heated(self):
+    return self._heated
+ 
+  @property
+  def cooled(self):
+    return self._cooled
+
+  @property  
+  def floor_area(self):
+    if self._floor_area is None:
+      self._floor_area = 0
+      for s in self._surfaces:
+        if s.type == 'Ground':
+          self._floor_area += s.area
+    return self._floor_area
+
+  @property
+  def bounded(self) -> List[ThermalBoundary]:
+    return self._bounded
+
+  @bounded.setter
+  def bounded(self, value):
+    self._bounded = value
+
+  @property
+  def surfaces(self):
+    return self._surfaces
+
+  @property
+  def additional_thermal_bridge_u_value(self):
+    return self._additional_thermal_bridge_u_value
+
+  @additional_thermal_bridge_u_value.setter
+  def additional_thermal_bridge_u_value(self, value):
+    self._additional_thermal_bridge_u_value = value
+
+  @property
+  def effective_thermal_capacity(self):
+    return self._effective_thermal_capacity
+
+  @effective_thermal_capacity.setter
+  def effective_thermal_capacity(self, value):
+    self._effective_thermal_capacity = value
+
+  @property
+  def indirectly_heated_area_ratio(self):
+    return self._indirectly_heated_area_ratio
+
+  @indirectly_heated_area_ratio.setter
+  def indirectly_heated_area_ratio(self, value):
+    self._indirectly_heated_area_ratio = value
+
+  @property
+  def infiltration_rate_system_on(self):
+    return self._infiltration_rate_system_on
+
+  @infiltration_rate_system_on.setter
+  def infiltration_rate_system_on(self, value):
+    self._infiltration_rate_system_on = value
+
+  @property
+  def infiltration_rate_system_off(self):
+    return self._infiltration_rate_system_off
+
+  @infiltration_rate_system_off.setter
+  def infiltration_rate_system_off(self, value):
+    self._infiltration_rate_system_off = value
+
diff --git a/city_model_structure/usage_zone.py b/city_model_structure/usage_zone.py
new file mode 100644
index 00000000..815fc407
--- /dev/null
+++ b/city_model_structure/usage_zone.py
@@ -0,0 +1,78 @@
+from city_model_structure.internal_gains import InternalGains
+from typing import List
+
+
+class UsageZone:
+  def __init__(self):
+    self._usage = None
+    self._internal_gains = None
+    self._heating_setpoint = None
+    self._heating_setback = None
+    self._cooling_setpoint = None
+    self._hours_day = None
+    self._days_year = None
+    self._mechanical_air_change = None
+
+  @property
+  def internal_gains(self) -> List[InternalGains]:
+    return self._internal_gains
+
+  @internal_gains.setter
+  def internal_gains(self, value):
+    self._internal_gains = value
+
+  @property
+  def heating_setpoint(self):
+    return self._heating_setpoint
+
+  @heating_setpoint.setter
+  def heating_setpoint(self, value):
+    self._heating_setpoint = value
+
+  @property
+  def heating_setback(self):
+    return self._heating_setback
+
+  @heating_setback.setter
+  def heating_setback(self, value):
+    self._heating_setback = value
+
+  @property
+  def cooling_setpoint(self):
+    return self._cooling_setpoint
+
+  @cooling_setpoint.setter
+  def cooling_setpoint(self, value):
+    self._cooling_setpoint = value
+
+  @property
+  def hours_day(self):
+    return self._hours_day
+
+  @hours_day.setter
+  def hours_day(self, value):
+    self._hours_day = value
+
+  @property
+  def days_year(self):
+    return self._days_year
+
+  @days_year.setter
+  def days_year(self, value):
+    self._days_year = value
+
+  @property
+  def mechanical_air_change(self):
+    return self._mechanical_air_change
+
+  @mechanical_air_change.setter
+  def mechanical_air_change(self, value):
+    self._mechanical_air_change = value
+
+  @property
+  def usage(self):
+    return self._usage
+
+  @usage.setter
+  def usage(self, value):
+    self._usage = value
diff --git a/city_model_structure/window.py b/city_model_structure/window.py
new file mode 100644
index 00000000..271a1aa9
--- /dev/null
+++ b/city_model_structure/window.py
@@ -0,0 +1,70 @@
+class Window:
+  def __init__(self):
+    # ToDo: construct this class
+    self._conductivity_wm_k = None
+    self._solar_transmittance_at_normal_incidence = None
+    self._front_side_solar_reflectance_at_normal_incidence = None
+    self._back_side_solar_reflectance_at_normal_incidence = None
+    self._frame_ratio = None
+    self._thickness_m = None
+    self._shgc = None
+
+  @property
+  def conductivity_wm_k(self):
+    return self._conductivity_wm_k
+
+  @conductivity_wm_k.setter
+  def conductivity_wm_k(self, value):
+    self._conductivity_wm_k = value
+
+  @property
+  def solar_transmittance_at_normal_incidence(self):
+    return self._solar_transmittance_at_normal_incidence
+
+  @solar_transmittance_at_normal_incidence.setter
+  def solar_transmittance_at_normal_incidence(self, value):
+    self._solar_transmittance_at_normal_incidence = value
+
+  @property
+  def front_side_solar_reflectance_at_normal_incidence(self):
+    return self._front_side_solar_reflectance_at_normal_incidence
+
+  @front_side_solar_reflectance_at_normal_incidence.setter
+  def front_side_solar_reflectance_at_normal_incidence(self, value):
+    self._front_side_solar_reflectance_at_normal_incidence = value
+
+  @property
+  def back_side_solar_reflectance_at_normal_incidence(self):
+    return self._back_side_solar_reflectance_at_normal_incidence
+
+  @back_side_solar_reflectance_at_normal_incidence.setter
+  def back_side_solar_reflectance_at_normal_incidence(self, value):
+    self._back_side_solar_reflectance_at_normal_incidence = value
+
+  @property
+  def frame_ratio(self):
+    return self._frame_ratio
+
+  @frame_ratio.setter
+  def frame_ratio(self, value):
+    self._frame_ratio = value
+
+  @frame_ratio.setter
+  def frame_ratio(self, value):
+    self._frame_ratio = value
+
+  @property
+  def thickness_m(self):
+    return self._thickness_m
+
+  @thickness_m.setter
+  def thickness_m(self, value):
+    self._thickness_m = value
+
+  @property
+  def shgc(self):
+    return self._shgc
+
+  @shgc.setter
+  def shgc(self, value):
+    self._shgc = value
diff --git a/geometry/geometry_factory.py b/geometry/geometry_factory.py
new file mode 100644
index 00000000..a6de2c52
--- /dev/null
+++ b/geometry/geometry_factory.py
@@ -0,0 +1,24 @@
+from geometry.geometry_feeders.city_gml import CityGml
+from city_model_structure.city import City
+
+
+class GeometryFactory:
+  def __init__(self, file_type, path):
+    self._file_type = file_type.lower()
+    self._path = path
+
+  @property
+  def citygml(self):
+    return CityGml(self._path).city
+
+  @property
+  def geojson(self):
+    raise Exception('Not implemented')
+
+  @property
+  def bim(self):
+    raise Exception('Not implemented')
+
+  @property
+  def city(self) -> City:
+    return getattr(self, self._file_type, lambda: None)
diff --git a/geometry/geometry_feeders/city_gml.py b/geometry/geometry_feeders/city_gml.py
new file mode 100644
index 00000000..00e78a1a
--- /dev/null
+++ b/geometry/geometry_feeders/city_gml.py
@@ -0,0 +1,106 @@
+import xmltodict
+import numpy as np
+from city_model_structure.city import City
+from city_model_structure.city_object import CityObject
+from city_model_structure.surface import Surface
+from helpers.geometry import Geometry
+
+
+class CityGml:
+  def __init__(self, path):
+    self._city = None
+    with open(path) as gml:
+      # Clean the namespaces is an important task to prevent wrong ns:field due poor citygml implementations
+      self._gml = xmltodict.parse(gml.read(), process_namespaces=True, xml_attribs=True, namespaces={
+        'http://www.opengis.net/gml': None,
+        'http://www.w3.org/2001/XMLSchema-instance': None,
+        'urn:oasis:names:tc:ciq:xsdschema:xAL:2.0': None,
+        'http://www.w3.org/1999/xlink': None,
+        'http://www.opengis.net/citygml/relief/2.0': None,
+        'http://www.opengis.net/citygml/building/2.0': None,
+        'http://www.opengis.net/citygml/building/2.0 http://schemas.opengis.net/citygml/building/2.0/building.xsd '
+        'http://www.opengis.net/citygml/relief/2.0 http://schemas.opengis.net/citygml/relief/2.0/relief.xsd" '
+        'xmlns="http://www.opengis.net/citygml/2.0': None,
+        'http://www.opengis.net/citygml/2.0': None
+      }, force_list=('cityObjectMember', 'curveMember'))
+      self._cityObjects = None
+      self._geometry = Geometry()
+      envelope = self._gml['CityModel']['boundedBy']['Envelope']
+      if '#text' in envelope['lowerCorner']:
+        self._lower_corner = np.fromstring(envelope['lowerCorner']['#text'], dtype=float, sep=' ')
+        self._upper_corner = np.fromstring(envelope['upperCorner']['#text'], dtype=float, sep=' ')
+      else:
+        self._lower_corner = np.fromstring(envelope['lowerCorner'], dtype=float, sep=' ')
+        self._upper_corner = np.fromstring(envelope['upperCorner'], dtype=float, sep=' ')
+
+      self._srs_name = envelope['@srsName']
+
+  @property
+  def content(self):
+    return self._gml
+
+  @property
+  def city(self):
+    if self._city is None:
+      city_objects = []
+      for o in self._gml['CityModel']['cityObjectMember']:
+        lod = 0
+        if 'lod1Solid' in o['Building']:
+          lod += 1
+        if 'lod2Solid' in o['Building']:
+          lod += 2
+        if 'lod3Solid' in o['Building']:
+          lod += 4
+        if 'lod4Solid' in o['Building']:
+          lod += 8
+          # ToDo: this is specific for Lod1 need to be modeled for higher lod's and lod combinations
+        name = o['Building']['@id']
+        lod_str = 'lod' + str(lod) + 'Solid'
+        lod_terrain_str = 'lod' + str(lod) + 'TerrainIntersection'
+        try:
+          surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList']['#text'])
+                      for s in o['Building'][lod_str]['Solid']['exterior']['CompositeSurface']['surfaceMember']]
+        except TypeError:
+          surfaces = [Surface(s['Polygon']['exterior']['LinearRing']['posList'])
+                      for s in o['Building'][lod_str]['Solid']['exterior']['CompositeSurface']['surfaceMember']]
+
+        if lod_terrain_str in o['Building']:
+          try:
+            curves = [c['LineString']['posList']['#text']
+                      for c in o['Building'][lod_terrain_str]['MultiCurve']['curveMember']]
+          except TypeError:
+            curves = [c['LineString']['posList']
+                      for c in o['Building'][lod_terrain_str]['MultiCurve']['curveMember']]
+
+          terrains = []
+          for curve in curves:
+            curve_points = np.fromstring(curve, dtype=float, sep=' ')
+            curve_points = self._geometry.to_points_matrix(curve_points, True)
+            terrains.append(curve_points)
+        else:
+
+          terrains = []
+          for s in surfaces:
+            ground = True
+            ground_points = []
+            for row in s.points:
+              # all surface vertex are in the lower z
+              ground_point = [row[0], row[1], self._lower_corner[2]]
+              ground_points = np.concatenate((ground_points, ground_point), axis=None)
+              ground = ground and self._geometry.almost_equal(row, ground_point)
+            if ground:
+              ground_points = self._geometry.to_points_matrix(ground_points)
+              # Do i need to remove the duplicated?
+              terrains.append(ground_points)
+        year_of_construction = None
+        function = None
+        if 'yearOfConstruction' in o['Building']:
+          year_of_construction = o['Building']['yearOfConstruction']
+        if 'function' in o['Building']:
+          function = o['Building']['function']
+
+        city_objects.append(CityObject(name, lod, surfaces, terrains, year_of_construction, function))
+      self._city = City(self._lower_corner, self._upper_corner, city_objects, self._srs_name)
+    return self._city
+
+
diff --git a/physics/physics_factory.py b/physics/physics_factory.py
new file mode 100644
index 00000000..7ac3a327
--- /dev/null
+++ b/physics/physics_factory.py
@@ -0,0 +1,27 @@
+from physics.physics_feeders.us_new_york_city_physics_parameters import UsNewYorkCityPhysicsParameters
+from physics.physics_feeders.us_physics_parameters import UsPhysicsParameters
+
+
+class PhysicsFactory:
+  def __init__(self, handler, city):
+    self._handler = handler.lower().replace(' ', '_')
+    self._city = city
+    self.factory()
+
+  def us_new_york_city(self):
+    UsNewYorkCityPhysicsParameters(self._city)
+
+  def us(self):
+    UsPhysicsParameters(self._city)
+
+  def ca(self):
+    raise Exception('Not implemented')
+
+  def de(self):
+    raise Exception('Not implemented')
+
+  def es(self):
+    raise Exception('Not implemented')
+
+  def factory(self):
+    getattr(self, self._handler, lambda: None)()
diff --git a/physics/physics_feeders/helpers/pluto_to_function.py b/physics/physics_feeders/helpers/pluto_to_function.py
new file mode 100644
index 00000000..b97750e4
--- /dev/null
+++ b/physics/physics_feeders/helpers/pluto_to_function.py
@@ -0,0 +1,228 @@
+class PlutoToFunction:
+
+  building_function = {
+    'A0': 'single family house',
+    'A1': 'single family house',
+    'A2': 'single family house',
+    'A3': 'single family house',
+    'A4': 'single family house',
+    'A5': 'single family house',
+    'A6': 'single family house',
+    'A7': 'single family house',
+    'A8': 'single family house',
+    'A9': 'single family house',
+    'B1': 'multifamily house',
+    'B2': 'multifamily house',
+    'B3': 'multifamily house',
+    'B9': 'multifamily house',
+    'C0': 'residential',
+    'C1': 'residential',
+    'C2': 'residential',
+    'C3': 'residential',
+    'C4': 'residential',
+    'C5': 'residential',
+    'C6': 'residential',
+    'C7': 'residential',
+    'C8': 'residential',
+    'C9': 'residential',
+    'D0': 'residential',
+    'D1': 'residential',
+    'D2': 'residential',
+    'D3': 'residential',
+    'D4': 'residential',
+    'D5': 'residential',
+    'D6': 'residential',
+    'D7': 'residential',
+    'D8': 'residential',
+    'D9': 'residential',
+    'E1': 'warehouse',
+    'E3': 'warehouse',
+    'E4': 'warehouse',
+    'E5': 'warehouse',
+    'E7': 'warehouse',
+    'E9': 'warehouse',
+    'F1': 'warehouse',
+    'F2': 'warehouse',
+    'F4': 'warehouse',
+    'F5': 'warehouse',
+    'F8': 'warehouse',
+    'F9': 'warehouse',
+    'G0': 'office',
+    'G1': 'office',
+    'G2': 'office',
+    'G3': 'office',
+    'G4': 'office',
+    'G5': 'office',
+    'G6': 'office',
+    'G7': 'office',
+    'G8': 'office',
+    'G9': 'office',
+    'H1': 'hotel',
+    'H2': 'hotel',
+    'H3': 'hotel',
+    'H4': 'hotel',
+    'H5': 'hotel',
+    'H6': 'hotel',
+    'H7': 'hotel',
+    'H8': 'hotel',
+    'H9': 'hotel',
+    'HB': 'hotel',
+    'HH': 'hotel',
+    'HR': 'hotel',
+    'HS': 'hotel',
+    'I1': 'hospital',
+    'I2': 'outpatient',
+    'I3': 'outpatient',
+    'I4': 'residential',
+    'I5': 'outpatient',
+    'I6': 'outpatient',
+    'I7': 'outpatient',
+    'I9': 'outpatient',
+    'J1': 'large office',
+    'J2': 'large office',
+    'J3': 'large office',
+    'J4': 'large office',
+    'J5': 'large office',
+    'J6': 'large office',
+    'J7': 'large office',
+    'J8': 'large office',
+    'J9': 'large office',
+    'K1': 'strip mall',
+    'K2': 'strip mall',
+    'K3': 'strip mall',
+    'K4': 'residential',
+    'K5': 'restaurant',
+    'K6': 'commercial',
+    'K7': 'commercial',
+    'K8': 'commercial',
+    'K9': 'commercial',
+    'L1': 'residential',
+    'L2': 'residential',
+    'L3': 'residential',
+    'L8': 'residential',
+    'L9': 'residential',
+    'M1': 'large office',
+    'M2': 'large office',
+    'M3': 'large office',
+    'M4': 'large office',
+    'M9': 'large office',
+    'N1': 'residential',
+    'N2': 'residential',
+    'N3': 'residential',
+    'N4': 'residential',
+    'N9': 'residential',
+    'O1': 'office',
+    'O2': 'office',
+    'O3': 'office',
+    'O4': 'office',
+    'O5': 'office',
+    'O6': 'office',
+    'O7': 'office',
+    'O8': 'office',
+    'O9': 'office',
+    'P1': 'large office',
+    'P2': 'hotel',
+    'P3': 'office',
+    'P4': 'office',
+    'P5': 'office',
+    'P6': 'office',
+    'P7': 'large office',
+    'P8': 'large office',
+    'P9': 'office',
+    'Q0': 'office',
+    'Q1': 'office',
+    'Q2': 'office',
+    'Q3': 'office',
+    'Q4': 'office',
+    'Q5': 'office',
+    'Q6': 'office',
+    'Q7': 'office',
+    'Q8': 'office',
+    'Q9': 'office',
+    'R0': 'residential',
+    'R1': 'residential',
+    'R2': 'residential',
+    'R3': 'residential',
+    'R4': 'residential',
+    'R5': 'residential',
+    'R6': 'residential',
+    'R7': 'residential',
+    'R8': 'residential',
+    'R9': 'residential',
+    'RA': 'residential',
+    'RB': 'residential',
+    'RC': 'residential',
+    'RD': 'residential',
+    'RG': 'residential',
+    'RH': 'residential',
+    'RI': 'residential',
+    'RK': 'residential',
+    'RM': 'residential',
+    'RR': 'residential',
+    'RS': 'residential',
+    'RW': 'residential',
+    'RX': 'residential',
+    'RZ': 'residential',
+    'S0': 'residential',
+    'S1': 'residential',
+    'S2': 'residential',
+    'S3': 'residential',
+    'S4': 'residential',
+    'S5': 'residential',
+    'S9': 'residential',
+    'T1': 'na',
+    'T2': 'na',
+    'T9': 'na',
+    'U0': 'warehouse',
+    'U1': 'warehouse',
+    'U2': 'warehouse',
+    'U3': 'warehouse',
+    'U4': 'warehouse',
+    'U5': 'warehouse',
+    'U6': 'warehouse',
+    'U7': 'warehouse',
+    'U8': 'warehouse',
+    'U9': 'warehouse',
+    'V0': 'na',
+    'V1': 'na',
+    'V2': 'na',
+    'V3': 'na',
+    'V4': 'na',
+    'V5': 'na',
+    'V6': 'na',
+    'V7': 'na',
+    'V8': 'na',
+    'V9': 'na',
+    'W1': 'primary school',
+    'W2': 'primary school',
+    'W3': 'secondary school',
+    'W4': 'secondary school',
+    'W5': 'secondary school',
+    'W6': 'secondary school',
+    'W7': 'secondary school',
+    'W8': 'primary school',
+    'W9': 'secondary school',
+    'Y1': 'large office',
+    'Y2': 'large office',
+    'Y3': 'large office',
+    'Y4': 'large office',
+    'Y5': 'large office',
+    'Y6': 'large office',
+    'Y7': 'large office',
+    'Y8': 'large office',
+    'Y9': 'large office',
+    'Z0': 'na',
+    'Z1': 'large office',
+    'Z2': 'na',
+    'Z3': 'na',
+    'Z4': 'na',
+    'Z5': 'na',
+    'Z6': 'na',
+    'Z7': 'na',
+    'Z8': 'na',
+    'Z9': 'na'
+  }
+  
+  @staticmethod
+  def function(pluto):
+    return PlutoToFunction.building_function[pluto]
diff --git a/physics/physics_feeders/helpers/us_to_library_types.py b/physics/physics_feeders/helpers/us_to_library_types.py
new file mode 100644
index 00000000..c7ed1e46
--- /dev/null
+++ b/physics/physics_feeders/helpers/us_to_library_types.py
@@ -0,0 +1,57 @@
+class UsToLibraryTypes(object):
+  standards = {
+    'ASHRAE Std189': 1,
+    'ASHRAE 90.1-2004': 2
+  }
+
+  reference_city_climate_zone = {
+    'Miami': 'ASHRAE_2004:1A',
+    'Houston': 'ASHRAE_2004:2A',
+    'Phoenix': 'ASHRAE_2004:2B',
+    'Atlanta': 'ASHRAE_2004:3A',
+    'Los Angeles': 'ASHRAE_2004:3B',
+    'Las Vegas': 'ASHRAE_2004:3B',
+    'San Francisco': 'ASHRAE_2004:3C',
+    'Baltimore': 'ASHRAE_2004:4A',
+    'Albuquerque': 'ASHRAE_2004:4B',
+    'Seattle': 'ASHRAE_2004:4C',
+    'Chicago': 'ASHRAE_2004:5A',
+    'Boulder': 'ASHRAE_2004:5B',
+    'Minneapolis': 'ASHRAE_2004:6A',
+    'Helena': 'ASHRAE_2004:6B',
+    'Duluth': 'ASHRAE_2004:7A',
+    'Fairbanks': 'ASHRAE_2004:8A'
+  }
+
+  window_types = ['window', 'door', 'skylight']
+
+  construction_types = {
+    'Wall': 'exterior wall',
+    'interior wall': 'interior wall',
+    'ground wall': 'ground wall',
+    'Ground': 'exterior slab',
+    'attic floor': 'attic floor',
+    'interior slab': 'interior slab',
+    'Roof': 'roof'
+  }
+
+  @staticmethod
+  def yoc_to_standard(year_of_construction):
+    if int(year_of_construction) < 2009:
+      standard = 'ASHRAE 90.1_2004'
+    else:
+      standard = 'ASHRAE 189.1_2009'
+    return standard
+
+  @staticmethod
+  def city_to_reference_city(city):
+    # ToDo: Dummy function that need to be implemented
+    reference_city = 'Baltimore'
+    if city is not None:
+      reference_city = 'Baltimore'
+    return reference_city
+
+  @staticmethod
+  def city_to_climate_zone(city):
+    reference_city = UsToLibraryTypes.city_to_reference_city(city)
+    return UsToLibraryTypes.reference_city_climate_zone[reference_city]
diff --git a/physics/physics_feeders/us_base_physics_parameters.py b/physics/physics_feeders/us_base_physics_parameters.py
new file mode 100644
index 00000000..bfd00443
--- /dev/null
+++ b/physics/physics_feeders/us_base_physics_parameters.py
@@ -0,0 +1,109 @@
+import xmltodict
+from city_model_structure.layer import Layer
+from city_model_structure.material import Material
+from pathlib import Path
+from physics.physics_feeders.helpers.us_to_library_types import UsToLibraryTypes
+
+
+class UsBasePhysicsParameters:
+  def __init__(self, climate_zone, city_objects, function_to_type):
+    self._climate_zone = climate_zone
+    self._city_objects = city_objects
+
+    # load US Library
+    path = str(Path.cwd() / 'data/physics/us_constructions.xml')
+    with open(path) as xml:
+      self._library = xmltodict.parse(xml.read(), force_list='layer')
+
+    # load US Archetypes
+    path = str(Path.cwd() / 'data/physics/us_archetypes.xml')
+    with open(path) as xml:
+      self._archetypes = xmltodict.parse(xml.read(), force_list='layer')
+    for city_object in self._city_objects:
+      building_type = function_to_type(city_object.function)
+      if building_type is None:
+        return
+      archetype = self.search_archetype(building_type,
+                                        UsToLibraryTypes.yoc_to_standard(city_object.year_of_construction),
+                                        self._climate_zone)
+      # ToDo:remove this in the future
+      # ToDo: Raise WrongArchetype if not all the surface types are defined for the given city_object
+      if archetype is None:
+        print(building_type, UsToLibraryTypes.yoc_to_standard(city_object.year_of_construction),
+              self._climate_zone)
+        raise Exception('Archetype not found for building')
+
+      city_object.average_storey_height = archetype['average_storey_height']['#text']
+      city_object.storeys_above_ground = archetype['number_of_storeys']['#text']
+
+      for thermal_zone in city_object.thermal_zones:
+        thermal_zone.effective_thermal_capacity = archetype['thermal_capacity']['#text']
+        thermal_zone.additional_thermal_bridge_u_value = archetype['extra_loses_due_to_thermal_bridges']['#text']
+        thermal_zone.indirectly_heated_area_ratio = archetype['indirect_heated_ratio']['#text']
+        thermal_zone.infiltration_rate_system_off = archetype['infiltration_rate_for_ventilation_system_off']['#text']
+        thermal_zone.infiltration_rate_system_on = archetype['infiltration_rate_for_ventilation_system_on']['#text']
+        for thermal_boundary in thermal_zone.bounded:
+          construction_type = UsToLibraryTypes.construction_types[thermal_boundary.type]
+          construction = UsBasePhysicsParameters.search_construction_in_archetype(archetype, construction_type)
+          construction_id = construction['@id']
+          c_lib = self.search_construction_type('construction', construction_id)
+          if 'outside_solar_absorptance' in c_lib:
+            thermal_boundary.outside_solar_absorptance = c_lib['outside_solar_absorptance']['#text']
+            thermal_boundary.outside_thermal_absorptance = c_lib['outside_thermal_absorptance']['#text']
+            thermal_boundary.outside_visible_absorptance = c_lib['outside_visible_absorptance']['#text']
+          thermal_boundary.window_ratio = construction['window_ratio']['#text']
+          thermal_boundary.layers = []
+          for current_layer in c_lib['layers']['layer']:
+            layer = Layer()
+            if 'thickness' in current_layer:
+              layer.thickness_m = current_layer['thickness']['#text']
+            material_lib = self.search_construction_type('material', current_layer['material'])
+            material = Material()
+            if 'conductivity' in material_lib:
+              material.conductivity_wm_k = material_lib['conductivity']['#text']
+              material.specific_heat_jkg_k = material_lib['specific_heat']['#text']
+              material.density_kg_m3 = material_lib['density']['#text']
+            material.solar_absorptance = material_lib['solar_absorptance']['#text']
+            material.thermal_absorptance = material_lib['thermal_absorptance']['#text']
+            material.visible_absorptance = material_lib['visible_absorptance']['#text']
+            material.no_mass = 'no_mass' in material_lib
+            if 'thermal_resistance' in material_lib:
+              material.thermal_resistance_m2k_w = material_lib['thermal_resistance']['#text']
+            else:
+              material.thermal_resistance_m2k_w = None
+            layer.material = material
+            thermal_boundary.layers.append(layer)
+          for opening in thermal_boundary.thermal_openings:
+            if construction['window'] is None:
+              continue
+            w_lib = self.search_construction_type('window', construction['window'])
+            opening.area = '0.0'
+            opening.conductivity_w_mk = w_lib['conductivity']['#text']
+            opening.frame_ratio = w_lib['frame_ratio']['#text']
+            opening.g_value = w_lib['solar_transmittance_at_normal_incidence']['#text']
+            opening.thickness_m = w_lib['thickness']['#text']
+            opening.inside_reflectance = 1.0-float(w_lib['back_side_solar_transmittance_at_normal_incidence']['#text'])
+            opening.outside_reflectance = \
+                1.0-float(w_lib['front_side_solar_transmittance_at_normal_incidence']['#text'])
+
+  def search_archetype(self, building_type, standard, climate_zone):
+    for archetype in self._archetypes['archetypes']['archetype']:
+      a_yc = str(archetype['@reference_standard'])
+      a_bt = str(archetype['@building_type'])
+      a_cz = str(archetype['@climate_zone'])
+      if (a_yc == str(standard)) and (a_bt == str(building_type)) and (a_cz == str(climate_zone)):
+        return archetype
+    return None
+
+  def search_construction_type(self, construction_type, construction_id):
+    for c_lib in self._library['library'][construction_type + 's'][construction_type]:
+      if construction_id == c_lib['@id']:
+        return c_lib
+    raise Exception('Archetype definition contains elements that does not exist in the library')
+
+  @staticmethod
+  def search_construction_in_archetype(archetype, construction_type):
+    for construction in archetype['constructions']['construction']:
+      if construction['@type'] == construction_type:
+        return construction
+    raise Exception('Construction type not found')
diff --git a/physics/physics_feeders/us_new_york_city_physics_parameters.py b/physics/physics_feeders/us_new_york_city_physics_parameters.py
new file mode 100644
index 00000000..c253b5ab
--- /dev/null
+++ b/physics/physics_feeders/us_new_york_city_physics_parameters.py
@@ -0,0 +1,10 @@
+from physics.physics_feeders.us_base_physics_parameters import UsBasePhysicsParameters
+from physics.physics_feeders.helpers.pluto_to_function import PlutoToFunction as Pf
+
+
+class UsNewYorkCityPhysicsParameters(UsBasePhysicsParameters):
+  def __init__(self, city):
+    self._city = city
+    climate_zone = 'ASHRAE_2004:4A'
+    super().__init__(climate_zone, self._city.city_objects, Pf.function)
+
diff --git a/physics/physics_feeders/us_physics_parameters.py b/physics/physics_feeders/us_physics_parameters.py
new file mode 100644
index 00000000..1a3c23e7
--- /dev/null
+++ b/physics/physics_feeders/us_physics_parameters.py
@@ -0,0 +1,11 @@
+from physics.physics_feeders.us_base_physics_parameters import UsBasePhysicsParameters
+from physics.physics_feeders.helpers.pluto_to_function import PlutoToFunction as Pf
+from physics.physics_feeders.helpers.us_to_library_types import UsToLibraryTypes
+
+
+class UsPhysicsParameters(UsBasePhysicsParameters):
+  def __init__(self, city):
+    self._city = city
+    self._climate_zone = UsToLibraryTypes.city_to_climate_zone(city.city_name)
+    super().__init__(self._climate_zone, self._city.city_objects, Pf.building_type_to_nrel)
+
diff --git a/usage/usage_factory.py b/usage/usage_factory.py
new file mode 100644
index 00000000..e287156f
--- /dev/null
+++ b/usage/usage_factory.py
@@ -0,0 +1,24 @@
+from usage.usage_feeders.de_usage_parameters import DeUsageParameters
+from usage.usage_feeders.us_new_york_city_usage_parameters import UsNewYorkCityUsageParameters
+
+
+class UsageFactory:
+  def __init__(self, handler, city):
+    self._handler = handler.lower().replace(' ', '_')
+    self._city = city
+    self.factory()
+
+  def us_new_york_city(self):
+    UsNewYorkCityUsageParameters(self._city)
+
+  def ca(self):
+    raise Exception('Not implemented')
+
+  def de(self):
+    DeUsageParameters(self._city)
+
+  def es(self):
+    raise Exception('Not implemented')
+
+  def factory(self):
+    getattr(self, self._handler, lambda: None)()
diff --git a/usage/usage_feeders/de_usage_parameters.py b/usage/usage_feeders/de_usage_parameters.py
new file mode 100644
index 00000000..559bfde4
--- /dev/null
+++ b/usage/usage_feeders/de_usage_parameters.py
@@ -0,0 +1,15 @@
+import xmltodict
+from pathlib import Path
+from usage.usage_feeders.helpers.us_function_to_usage import UsFunctionToUsage
+
+
+class DeUsageParameters:
+  def __init__(self, city_objects):
+    self._city_objects = city_objects
+
+    # load US Library
+    path = str(Path.cwd() / 'data/usage/de_library.xml')
+    with open(path) as xml:
+      self._library = xmltodict.parse(xml.read())
+    for city_object in city_objects:
+       UsFunctionToUsage.function_to_usage(city_object.function)
diff --git a/usage/usage_feeders/helpers/us_function_to_usage.py b/usage/usage_feeders/helpers/us_function_to_usage.py
new file mode 100644
index 00000000..fc3f3019
--- /dev/null
+++ b/usage/usage_feeders/helpers/us_function_to_usage.py
@@ -0,0 +1,251 @@
+class UsFunctionToUsage:
+  building_usage = {
+    'full service restaurant': 'restaurant',
+    'highrise apartment': 'residential',
+    'hospital': 'health care',
+    'large hotel': 'hotel',
+    'large office': 'office and administration',
+    'medium office': 'office and administration',
+    'midrise apartment': 'residential',
+    'outpatient healthcare': 'health care',
+    'primary school': 'education',
+    'quick service restaurant': 'restaurant',
+    'secondary school': 'education',
+    'small hotel': 'hotel',
+    'small office': 'office and administration',
+    'stand alone retail': 'retail',
+    'strip mall': 'hall',
+    'supermarket': 'retail',
+    'warehouse': 'industry'
+  }
+
+  building_function = {
+    'A0': 'single family house',
+    'A1': 'single family house',
+    'A2': 'single family house',
+    'A3': 'single family house',
+    'A4': 'single family house',
+    'A5': 'single family house',
+    'A6': 'single family house',
+    'A7': 'single family house',
+    'A8': 'single family house',
+    'A9': 'single family house',
+    'B1': 'multifamily house',
+    'B2': 'multifamily house',
+    'B3': 'multifamily house',
+    'B9': 'multifamily house',
+    'C0': 'residential',
+    'C1': 'residential',
+    'C2': 'residential',
+    'C3': 'residential',
+    'C4': 'residential',
+    'C5': 'residential',
+    'C6': 'residential',
+    'C7': 'residential',
+    'C8': 'residential',
+    'C9': 'residential',
+    'D0': 'residential',
+    'D1': 'residential',
+    'D2': 'residential',
+    'D3': 'residential',
+    'D4': 'residential',
+    'D5': 'residential',
+    'D6': 'residential',
+    'D7': 'residential',
+    'D8': 'residential',
+    'D9': 'residential',
+    'E1': 'warehouse',
+    'E3': 'warehouse',
+    'E4': 'warehouse',
+    'E5': 'warehouse',
+    'E7': 'warehouse',
+    'E9': 'warehouse',
+    'F1': 'warehouse',
+    'F2': 'warehouse',
+    'F4': 'warehouse',
+    'F5': 'warehouse',
+    'F8': 'warehouse',
+    'F9': 'warehouse',
+    'G0': 'office',
+    'G1': 'office',
+    'G2': 'office',
+    'G3': 'office',
+    'G4': 'office',
+    'G5': 'office',
+    'G6': 'office',
+    'G7': 'office',
+    'G8': 'office',
+    'G9': 'office',
+    'H1': 'hotel',
+    'H2': 'hotel',
+    'H3': 'hotel',
+    'H4': 'hotel',
+    'H5': 'hotel',
+    'H6': 'hotel',
+    'H7': 'hotel',
+    'H8': 'hotel',
+    'H9': 'hotel',
+    'HB': 'hotel',
+    'HH': 'hotel',
+    'HR': 'hotel',
+    'HS': 'hotel',
+    'I1': 'hospital',
+    'I2': 'outpatient',
+    'I3': 'outpatient',
+    'I4': 'residential',
+    'I5': 'outpatient',
+    'I6': 'outpatient',
+    'I7': 'outpatient',
+    'I9': 'outpatient',
+    'J1': 'large office',
+    'J2': 'large office',
+    'J3': 'large office',
+    'J4': 'large office',
+    'J5': 'large office',
+    'J6': 'large office',
+    'J7': 'large office',
+    'J8': 'large office',
+    'J9': 'large office',
+    'K1': 'strip mall',
+    'K2': 'strip mall',
+    'K3': 'strip mall',
+    'K4': 'residential',
+    'K5': 'restaurant',
+    'K6': 'commercial',
+    'K7': 'commercial',
+    'K8': 'commercial',
+    'K9': 'commercial',
+    'L1': 'residential',
+    'L2': 'residential',
+    'L3': 'residential',
+    'L8': 'residential',
+    'L9': 'residential',
+    'M1': 'large office',
+    'M2': 'large office',
+    'M3': 'large office',
+    'M4': 'large office',
+    'M9': 'large office',
+    'N1': 'residential',
+    'N2': 'residential',
+    'N3': 'residential',
+    'N4': 'residential',
+    'N9': 'residential',
+    'O1': 'office',
+    'O2': 'office',
+    'O3': 'office',
+    'O4': 'office',
+    'O5': 'office',
+    'O6': 'office',
+    'O7': 'office',
+    'O8': 'office',
+    'O9': 'office',
+    'P1': 'large office',
+    'P2': 'hotel',
+    'P3': 'office',
+    'P4': 'office',
+    'P5': 'office',
+    'P6': 'office',
+    'P7': 'large office',
+    'P8': 'large office',
+    'P9': 'office',
+    'Q0': 'office',
+    'Q1': 'office',
+    'Q2': 'office',
+    'Q3': 'office',
+    'Q4': 'office',
+    'Q5': 'office',
+    'Q6': 'office',
+    'Q7': 'office',
+    'Q8': 'office',
+    'Q9': 'office',
+    'R0': 'residential',
+    'R1': 'residential',
+    'R2': 'residential',
+    'R3': 'residential',
+    'R4': 'residential',
+    'R5': 'residential',
+    'R6': 'residential',
+    'R7': 'residential',
+    'R8': 'residential',
+    'R9': 'residential',
+    'RA': 'residential',
+    'RB': 'residential',
+    'RC': 'residential',
+    'RD': 'residential',
+    'RG': 'residential',
+    'RH': 'residential',
+    'RI': 'residential',
+    'RK': 'residential',
+    'RM': 'residential',
+    'RR': 'residential',
+    'RS': 'residential',
+    'RW': 'residential',
+    'RX': 'residential',
+    'RZ': 'residential',
+    'S0': 'residential',
+    'S1': 'residential',
+    'S2': 'residential',
+    'S3': 'residential',
+    'S4': 'residential',
+    'S5': 'residential',
+    'S9': 'residential',
+    'T1': 'na',
+    'T2': 'na',
+    'T9': 'na',
+    'U0': 'warehouse',
+    'U1': 'warehouse',
+    'U2': 'warehouse',
+    'U3': 'warehouse',
+    'U4': 'warehouse',
+    'U5': 'warehouse',
+    'U6': 'warehouse',
+    'U7': 'warehouse',
+    'U8': 'warehouse',
+    'U9': 'warehouse',
+    'V0': 'na',
+    'V1': 'na',
+    'V2': 'na',
+    'V3': 'na',
+    'V4': 'na',
+    'V5': 'na',
+    'V6': 'na',
+    'V7': 'na',
+    'V8': 'na',
+    'V9': 'na',
+    'W1': 'primary school',
+    'W2': 'primary school',
+    'W3': 'secondary school',
+    'W4': 'secondary school',
+    'W5': 'secondary school',
+    'W6': 'secondary school',
+    'W7': 'secondary school',
+    'W8': 'primary school',
+    'W9': 'secondary school',
+    'Y1': 'large office',
+    'Y2': 'large office',
+    'Y3': 'large office',
+    'Y4': 'large office',
+    'Y5': 'large office',
+    'Y6': 'large office',
+    'Y7': 'large office',
+    'Y8': 'large office',
+    'Y9': 'large office',
+    'Z0': 'na',
+    'Z1': 'large office',
+    'Z2': 'na',
+    'Z3': 'na',
+    'Z4': 'na',
+    'Z5': 'na',
+    'Z6': 'na',
+    'Z7': 'na',
+    'Z8': 'na',
+    'Z9': 'na'
+  }
+
+  @staticmethod
+  def function_to_usage(building_function):
+    return UsFunctionToUsage.building_usage[building_function]
+
+  @staticmethod
+  def pluto_to_function(pluto):
+    return UsFunctionToUsage.building_function[pluto]
diff --git a/usage/usage_feeders/us_base_usage_parameters.py b/usage/usage_feeders/us_base_usage_parameters.py
new file mode 100644
index 00000000..e982087b
--- /dev/null
+++ b/usage/usage_feeders/us_base_usage_parameters.py
@@ -0,0 +1,78 @@
+import xmltodict
+from pathlib import Path
+from city_model_structure.usage_zone import UsageZone
+from city_model_structure.internal_gains import InternalGains
+
+
+class UsBaseUsageParameters:
+  def __init__(self, city, function_to_usage):
+    self._city = city
+    # ToDo: this is using the german library as a temporary approach, need to use/define a library for US
+    path = str(Path.cwd() / 'data/usage/de_library.xml')
+    with open(path) as xml:
+      self._library = xmltodict.parse(xml.read(), force_list='zoneUsageVariant')
+    for city_object in self._city.city_objects:
+      #ToDo: Right now is just one usage zone but will be multiple in the future
+      usage_zone = UsageZone()
+      usage_zone.usage = function_to_usage(city_object.function)
+      for zone_usage_type in self._library['buildingUsageLibrary']['zoneUsageType']:
+        if zone_usage_type['id'] != usage_zone.usage:
+          if 'zoneUsageVariant' in zone_usage_type:
+            for usage_zone_variant in zone_usage_type['zoneUsageVariant']:
+              if usage_zone_variant['id'] == usage_zone.usage:
+                # pre-initialize the usage zone with the main type
+                usage_zone = UsBaseUsageParameters._parse_zone_usage_type(zone_usage_type, usage_zone)
+                usage_zone = UsBaseUsageParameters._parse_zone_usage_variant(usage_zone_variant, usage_zone)
+                city_object.usage_zone = [usage_zone]
+                break
+          continue
+        else:
+          city_object.usage_zone = [UsBaseUsageParameters._parse_zone_usage_type(zone_usage_type, usage_zone)]
+          break
+      if city_object.usage_zone is None:
+        print(city_object.function)
+        raise Exception('Usage not found for building function')
+
+
+  @staticmethod
+  def _parse_zone_usage_type(zone_usage_type, usage_zone):
+    usage_zone.hours_day = zone_usage_type['occupancy']['usageHoursPerDay']
+    usage_zone.days_year = zone_usage_type['occupancy']['usageDaysPerYear']
+    usage_zone.cooling_setpoint = zone_usage_type['endUses']['space_cooling']['coolingSetPointTemperature']
+    usage_zone.heating_setpoint = zone_usage_type['endUses']['space_heating']['heatingSetPointTemperature']
+    usage_zone.heating_setback = zone_usage_type['endUses']['space_heating']['heatingSetBackTemperature']
+    if 'ventilation' in zone_usage_type['endUses'] and zone_usage_type['endUses']['ventilation'] is not None:
+      usage_zone.mechanical_air_change = zone_usage_type['endUses']['ventilation']['mechanicalAirChangeRate']
+    int_gains = InternalGains()
+    int_gains.latent_fraction = zone_usage_type['occupancy']['internGains']['latentFraction']
+    int_gains.convective_fraction = zone_usage_type['occupancy']['internGains']['convectiveFraction']
+    int_gains.average_internal_gain_w_m2 = zone_usage_type['occupancy']['internGains']['averageInternGainPerSqm']
+    int_gains.radiative_fraction = zone_usage_type['occupancy']['internGains']['radiantFraction']
+    usage_zone.internal_gains = [int_gains]
+    return usage_zone
+
+  @staticmethod
+  def _parse_zone_usage_variant(usage_zone_variant, usage_zone):
+    # for the variants all is optional because it mimics the inheritance concept from OOP
+    if 'usageHoursPerDay' in usage_zone_variant['occupancy']:
+      usage_zone.hours_day = usage_zone_variant['occupancy']['usageHoursPerDay']
+    if 'usageDaysPerYear' in usage_zone_variant['occupancy']:
+      usage_zone.days_year = usage_zone_variant['occupancy']['usageDaysPerYear']
+    if 'coolingSetPointTemperature' in usage_zone_variant['endUses']['space_cooling']:
+      usage_zone.cooling_setpoint = usage_zone_variant['endUses']['space_cooling']['coolingSetPointTemperature']
+    if 'heatingSetPointTemperature' in usage_zone_variant['endUses']['space_heating']:
+      usage_zone.heating_setpoint = usage_zone_variant['endUses']['space_heating']['heatingSetPointTemperature']
+    if 'heatingSetBackTemperature' in usage_zone_variant['endUses']['space_heating']:
+      usage_zone.heating_setback = usage_zone_variant['endUses']['space_heating']['heatingSetBackTemperature']
+    if 'ventilation' in usage_zone_variant['endUses'] and usage_zone_variant['endUses']['ventilation'] is not None:
+      usage_zone.mechanical_air_change = usage_zone_variant['endUses']['ventilation']['mechanicalAirChangeRate']
+    if 'latentFraction' in usage_zone_variant['occupancy']['internGains']:
+      usage_zone.int_gains[0].latent_fraction = usage_zone_variant['occupancy']['internGains']['latentFraction']
+    if 'convectiveFraction' in usage_zone_variant['occupancy']['internGains']:
+      usage_zone.int_gains[0].convective_fraction = usage_zone_variant['occupancy']['internGains']['convectiveFraction']
+    if 'averageInternGainPerSqm' in usage_zone_variant['occupancy']['internGains']:
+      usage_zone.int_gains[0].average_internal_gain_w_m2 = \
+        usage_zone_variant['occupancy']['internGains']['averageInternGainPerSqm']
+    if 'radiantFraction' in usage_zone_variant['occupancy']['internGains']:
+      usage_zone.int_gains[0].radiative_fraction = usage_zone_variant['occupancy']['internGains']['radiantFraction']
+    return usage_zone
diff --git a/usage/usage_feeders/us_new_york_city_usage_parameters.py b/usage/usage_feeders/us_new_york_city_usage_parameters.py
new file mode 100644
index 00000000..dc5d5ad1
--- /dev/null
+++ b/usage/usage_feeders/us_new_york_city_usage_parameters.py
@@ -0,0 +1,7 @@
+from usage.usage_feeders.us_base_usage_parameters import UsBaseUsageParameters
+from usage.usage_feeders.helpers.us_function_to_usage import UsFunctionToUsage
+
+
+class UsNewYorkCityUsageParameters(UsBaseUsageParameters):
+  def __init__(self, city):
+    super().__init__(city, UsFunctionToUsage.pluto_to_function)
\ No newline at end of file
diff --git a/usage/usage_feeders/us_usage_parameters.py b/usage/usage_feeders/us_usage_parameters.py
new file mode 100644
index 00000000..bd645cd7
--- /dev/null
+++ b/usage/usage_feeders/us_usage_parameters.py
@@ -0,0 +1,7 @@
+from usage.usage_feeders.us_base_usage_parameters import UsBaseUsageParameters
+from usage.usage_feeders.helpers.us_function_to_usage import UsFunctionToUsage
+
+
+class UsUsageParameters(UsBaseUsageParameters):
+  def __init__(self, city):
+    super().__init__(city, UsFunctionToUsage.function_to_usage)
diff --git a/usage/usagelibrary.py b/usage/usagelibrary.py
new file mode 100644
index 00000000..fd0d3432
--- /dev/null
+++ b/usage/usagelibrary.py
@@ -0,0 +1,10 @@
+class UsageLibrary:
+  def __init__(self):
+    self.standard = 1
+    self.internal_gains = [4.2]
+    self.heating_setpoint = [20.0]
+    self.heating_setback = [16.0]
+    self.cooling_setpoint = [25.0]
+    self.hours_day = [17]
+    self.days_year = [365]
+    self.min_air_change = [0.34]