diff --git a/city_model_structure/city_object.py b/city_model_structure/city_object.py index 5d368ff1..5babaf5d 100644 --- a/city_model_structure/city_object.py +++ b/city_model_structure/city_object.py @@ -15,7 +15,7 @@ from shapely.geometry import MultiPolygon import numpy as np import matplotlib.patches as patches -from helpers.geometry import Geometry +from helpers.geometry_helper import GeometryHelper from city_model_structure.usage_zone import UsageZone @@ -35,7 +35,7 @@ class CityObject: self._year_of_construction = year_of_construction self._function = function self._lower_corner = lower_corner - self._geometry = Geometry() + self._geometry = GeometryHelper() self._average_storey_height = None self._storeys_above_ground = None self._foot_print = None diff --git a/city_model_structure/polyhedron.py b/city_model_structure/polyhedron.py index 9fdaa2cd..d1665029 100644 --- a/city_model_structure/polyhedron.py +++ b/city_model_structure/polyhedron.py @@ -5,7 +5,7 @@ Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@conc """ import numpy as np import stl -from helpers.geometry import Geometry +from helpers.geometry_helper import GeometryHelper class Polyhedron: @@ -20,7 +20,7 @@ class Polyhedron: self._faces = None self._vertices = None self._mesh = None - self._geometry = Geometry() + self._geometry = GeometryHelper() def _position_of(self, point): vertices = self.vertices diff --git a/config/configuration.ini b/config/configuration.ini index 9322833d..81cdb84f 100644 --- a/config/configuration.ini +++ b/config/configuration.ini @@ -14,6 +14,3 @@ indirectly_heated_area_ratio = 0 infiltration_rate_system_on = 0 outside_solar_absorptance = 0.2 shortwave_reflectance = 0.8 - - -#EOF \ No newline at end of file diff --git a/geometry/geometry_feeders/city_gml.py b/geometry/geometry_feeders/city_gml.py index b9948711..fead90b0 100644 --- a/geometry/geometry_feeders/city_gml.py +++ b/geometry/geometry_feeders/city_gml.py @@ -8,7 +8,7 @@ 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 +from helpers.geometry_helper import GeometryHelper class CityGml: @@ -32,7 +32,7 @@ class CityGml: 'http://www.opengis.net/citygml/2.0': None }, force_list=('cityObjectMember', 'curveMember')) self._cityObjects = None - self._geometry = Geometry() + self._geometry = GeometryHelper() envelope = self._gml['CityModel']['boundedBy']['Envelope'] if '#text' in envelope['lowerCorner']: self._lower_corner = np.fromstring(envelope['lowerCorner']['#text'], dtype=float, sep=' ') diff --git a/helpers/configuration_helper.py b/helpers/configuration_helper.py new file mode 100644 index 00000000..02d4cd66 --- /dev/null +++ b/helpers/configuration_helper.py @@ -0,0 +1,91 @@ +""" +Configuration helper +SPDX - License - Identifier: LGPL - 3.0 - or -later +Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca +""" +import configparser +from pathlib import Path + + +class ConfigurationHelper: + """ + Configuration class + """ + def __init__(self): + base_path = Path().resolve().parent + config_file = Path(base_path / 'libs/config/configuration.ini').resolve() + self._config = configparser.ConfigParser() + self._config.read(config_file) + + @property + def h_i(self): + """ + Configured internal convective coefficient in W/m2K + :return: float + """ + return self._config.getfloat('convective_fluxes', 'h_i') + + @property + def h_e(self): + """ + Configured external convective coefficient in W/m2K + :return: float + """ + return self._config.getfloat('convective_fluxes', 'h_e') + + @property + def frame_ratio(self): + """ + Configured frame ratio + :return: float + """ + return self._config.getfloat('windows', 'frame_ratio') + + @property + def heated(self): + """ + Configured heated flag + :return: Boolean + """ + return self._config.getboolean('thermal_zones', 'heated') + + @property + def cooled(self): + """ + Configured cooled flag + :return: Boolean + """ + return self._config.getboolean('thermal_zones', 'cooled') + + @property + def additional_thermal_bridge_u_value(self): + """ + Configured additional thermal bridge u value W/m2K + :return: + """ + return self._config.getfloat('thermal_zones', 'additional_thermal_bridge_u_value') + + @property + def indirectly_heated_area_ratio(self): + """ + Configured indirectly heated area ratio + :return: float + """ + + return self._config.getfloat('thermal_zones', 'indirectly_heated_area_ratio') + + @property + def infiltration_rate_system_on(self): + """ + Configured infiltration rate system on in air change per hour + :return: float + """ + return self._config.getfloat('thermal_zones', 'infiltration_rate_system_on') + + @property + def outside_solar_absorptance(self): + """ + Configured infiltration rate system off in air change per hour + :return: float + """ + return self._config.getfloat('thermal_zones', 'outside_solar_absorptance') diff --git a/helpers/geometry_helper.py b/helpers/geometry_helper.py new file mode 100644 index 00000000..a2558392 --- /dev/null +++ b/helpers/geometry_helper.py @@ -0,0 +1,165 @@ +""" +Geometry helper +SPDX - License - Identifier: LGPL - 3.0 - or -later +Copyright © 2020 Project Author Guille Gutierrez guillermo.gutierrezmorote@concordia.ca +""" +import math +import numpy as np +from trimesh import Trimesh +from trimesh import intersections +import open3d as o3d + + +class GeometryHelper: + """ + Geometry helper class + """ + def __init__(self, delta=0.5): + self._delta = delta + + def almost_equal(self, v1, v2): + """ + Compare two points and decides if they are almost equal (quadratic error under delta) + :param v1: [x,y,z] + :param v2: [x,y,z] + :return: Boolean + """ + delta = math.sqrt(pow((v1[0] - v2[0]), 2) + pow((v1[1] - v2[1]), 2) + pow((v1[2] - v2[2]), 2)) + return delta <= self._delta + + def is_almost_same_surface(self, s1, s2): + """ + Compare two surfaces and decides if they are almost equal (quadratic error under delta) + :param s1: Surface + :param s2: Surface + :return: Boolean + """ + # delta is grads an need to be converted into radians + delta = np.rad2deg(self._delta) + difference = (s1.inclination - s2.inclination) % math.pi + if abs(difference) > delta: + return False + # s1 and s2 are at least almost parallel surfaces + # calculate distance point to plane using all the vertex + # select surface1 value for the point (X,Y,Z) where two of the values are 0 + minimum_distance = self._delta + 1 + parametric = s2.polygon.get_parametric() + n2 = s2.normal + for point in s1.points: + distance = abs( + (point[0] * parametric[0]) + (point[1] * parametric[1]) + (point[2] * parametric[2]) + parametric[3]) + normal_module = math.sqrt(pow(n2[0], 2) + pow(n2[1], 2) + pow(n2[2], 2)) + + if normal_module == 0: + continue + distance = distance / normal_module + if distance < minimum_distance: + minimum_distance = distance + if minimum_distance <= self._delta: + break + if minimum_distance > self._delta or s1.intersect(s2) is None: + return False + else: + return True + + @staticmethod + def to_points_matrix(points, remove_last=False): + """ + Transform a point vector into a point matrix + :param points: [x, y, z, x, y, z ...] + :param remove_last: Boolean + :return: [[x,y,z],[x,y,z]...] + """ + rows = points.size // 3 + points = points.reshape(rows, 3) + if remove_last: + points = np.delete(points, rows - 1, 0) + return points + + @staticmethod + def _segment_list_to_point_cloud(segment_list): + point_list = np.asarray(segment_list[0]) + for segment in segment_list: + for new_point in segment: + found = False + for point in point_list: + same_point = np.allclose(new_point, point) + if same_point: + found = True + break + if not found: + point_list = np.concatenate((point_list, [new_point])) + return point_list + + @staticmethod + def _point_cloud_to_mesh(point_list, normal_list): + # Return a mesh composed only by triangles + pcd = o3d.geometry.PointCloud() + pcd.points = o3d.utility.Vector3dVector(point_list) + pcd.normals = o3d.utility.Vector3dVector(normal_list) + distances = pcd.compute_nearest_neighbor_distance() + avg_dist = np.mean(distances) + radius = 3 * avg_dist + bpa_mesh = o3d.geometry.TriangleMesh().create_from_point_cloud_ball_pivoting( + pcd, o3d.utility.DoubleVector([radius, radius * 2])) + mesh_result = Trimesh(vertices=np.asarray(bpa_mesh.vertices), faces=np.asarray(bpa_mesh.triangles)) + return mesh_result + + @staticmethod + def _merge_meshes(mesh1, mesh2): + v_1 = mesh1.vertices + f_1 = mesh1.faces + v_2 = mesh2.vertices + f_2 = mesh2.faces + length = len(v_1) + v_merge = np.concatenate((v_1, v_2)) + f_merge = np.asarray(f_1) + + for item in f_2: + point1 = item.item(0) + length + point2 = item.item(1) + length + point3 = item.item(2) + length + surface = np.asarray([point1, point2, point3]) + f_merge = np.concatenate((f_merge, [surface])) + + mesh_merge = Trimesh(vertices=v_merge, faces=f_merge) + + return mesh_merge + + @staticmethod + def divide_mesh_by_plane(mesh, normal_plane, point_plane): + """ + Divide a mesh by a plane + :param mesh: Trimesh + :param normal_plane: [x, y, z] + :param point_plane: [x, y, z] + :return: [Trimesh] + """ + # The first mesh returns the positive side of the plane and the second the negative side. + # If the plane does not divide the mesh (i.e. it does not touch it or it is coplanar with one or more faces), + # then it returns only the original mesh. + normal_plane_opp = [None] * len(normal_plane) + for i in range(0, len(normal_plane)): + normal_plane_opp[i] = - normal_plane[i] + + normal = [normal_plane, normal_plane_opp] + normal_opp = [normal_plane_opp, normal_plane] + mesh_final = [] + for i in range(0, 2): + mesh_1 = intersections.slice_mesh_plane(mesh, normal[i], point_plane) + mesh_1_segments = intersections.mesh_plane(mesh, normal[i], point_plane) + boo = mesh.difference(mesh_1, engine='blender') + print(boo) + quit() + if len(mesh_1_segments) <= 0 or len(mesh_1.faces) == len(mesh.faces): + mesh_final.append(mesh) + break + else: + points = GeometryHelper._segment_list_to_point_cloud(mesh_1_segments) + points_normals = [[None] * 3] * len(points) + for j in range(0, len(points_normals)): + points_normals[j] = normal_opp[i] + mesh_2 = GeometryHelper._point_cloud_to_mesh(points, points_normals) + mesh_final.append(GeometryHelper._merge_meshes(mesh_1, mesh_2)) + + return mesh_final diff --git a/tests/test_geometry_factory.py b/tests/test_geometry_factory.py index 7d763035..77390fbc 100644 --- a/tests/test_geometry_factory.py +++ b/tests/test_geometry_factory.py @@ -23,7 +23,7 @@ class TestGeometryFactory(TestCase): def _get_citygml(self): if self._city_gml is None: - file_path = (self._example_path / 'buildings.gml').resolve() + file_path = (self._example_path / '2050 bp_2buildings.gml').resolve() self._city_gml = GeometryFactory('citygml', file_path).city self.assertIsNotNone(self._city_gml, 'city is none') return self._city_gml