Skip to content

Modules

Vertex

Vertex and Vertex related functions.

A Vertex is a 3-coordinate class, that can be used to represent a Point3d or a Vector3d Includes utilities to go from/to numpy array and OpenStudio's Point3d

Vertex

Point3d and Vector3d.

Source code in geomeffibem/vertex.py 
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
class Vertex:
    """Point3d and Vector3d."""

    @staticmethod
    def from_numpy(arr):
        """Factory method to construct from a numpy array (or a list) of 3 coords."""
        if not isinstance(arr, np.ndarray):
            arr = np.array(arr)
        if arr.shape != (3,):
            raise ValueError(f"Expected a numpy array with a dimension (3, ) (a Vector3d), got {arr.shape}")
        return Vertex(arr[0], arr[1], arr[2])

    @staticmethod
    def from_Point3d(pt: openstudio.Point3d):
        """Factory method to construct from an openstudio.Point3d."""
        return Vertex(pt.x(), pt.y(), pt.z())

    def __init__(self, x, y, z):
        """Vertex constructor."""
        self.x = x
        self.y = y
        self.z = z
        self.surface = None

    def copy(self):
        """Make a copy of this Vertex."""
        return Vertex(self.x, self.y, self.z)

    def get_coords_on_plane(self, plane='xy') -> Tuple[float, float]:
        """Returns two coordinates on a given plane."""
        if plane == 'xy':
            return self.x, self.y
        elif plane == 'xz':
            return self.x, self.z
        elif plane == 'yz':
            return self.y, self.z

        raise ValueError("Expected plane to be 'xy', 'xz' or 'yz'.")

    def length(self) -> float:
        """Get the length of the vector."""
        return np.sqrt(np.sum(self.to_numpy() ** 2))

    def normalize(self) -> Vertex:
        """Normalize to a length of 1, returns a copy."""
        v = self.copy()
        v.setLength(1.0)
        return v

    def setLength(self, newLength: float) -> None:
        """Change length of vector, in place."""
        currentLength = self.length()
        if currentLength > 0:
            mult = newLength / currentLength
            self.x *= mult
            self.y *= mult
            self.z *= mult
        else:
            raise ValueError("Cannot normalize a vector of length 0")

    def dot(self, other) -> float:
        """Computes the dot / scalar / inner product product of two vectors.

        (a.b).
        """
        # return np.dot(v.to_numpy(), v2.to_numpy())
        return self.x * other.x + self.y * other.y + self.z * other.z

    def cross(self, other, normalize: bool = False) -> Vertex:
        """Computes the cross product (a x b), which is a vector perpendicular to both a and b."""
        v = Vertex(
            x=(self.y * other.z - self.z * other.y),
            y=(self.z * other.x - self.x * other.z),
            z=(self.x * other.y - self.y * other.x),
        )
        if normalize:
            return v.normalize()
        return v

    def outer_product(self, other) -> np.ndarray:
        """Compute the outer product of this by another vector."""
        return np.outer(self.to_numpy(), other.to_numpy())

    def __add__(self, other) -> Vertex:
        """Return a + b."""
        return Vertex(x=self.x + other.x, y=self.y + other.y, z=self.z + other.z)

    def __neg__(self):
        """Return obj negated (-obj)."""
        return Vertex(-self.x, -self.y, -self.z)

    def __sub__(self, other) -> Vertex:
        """Return a - b."""
        return self + -other
        # return Vertex(
        #     x=self.x - other.x,
        #     y=self.y - other.y,
        #     z=self.z - other.z
        # )

    def __mul__(self, other) -> Vertex:
        """Multiplies each coordinate by a scalar."""
        if not isinstance(other, int) and not isinstance(other, float):
            raise ValueError("Multiplication of a vertex by something else than a numeric is not supported")
        return Vertex(self.x * other, self.y * other, self.z * other)

    def __rmul__(self, other) -> Vertex:
        """Multiplies each coordinate by a scalar."""
        if not isinstance(other, int) and not isinstance(other, float):
            raise ValueError("Multiplication of a vertex by something else than a numeric is not supported")
        return self.__mul__(other)

    def __truediv__(self, other) -> Vertex:
        """Divides each coordinate by a scalar."""
        if not isinstance(other, int) and not isinstance(other, float):
            raise ValueError("Division of a vertex by something else than a numeric is not supported")
        return Vertex(self.x / other, self.y / other, self.z / other)

    def __floordiv__(self, other) -> Vertex:
        """Divides each coordinate by a scalar."""
        if not isinstance(other, int) and not isinstance(other, float):
            raise ValueError("Division of a vertex by something else than a numeric is not supported")
        return Vertex(self.x // other, self.y // other, self.z // other)

    def to_numpy(self) -> np.ndarray:
        """Export to a numpy array of 3 coordinates."""
        return np.array([self.x, self.y, self.z])

    def to_Point3d(self) -> openstudio.Point3d:
        """Export to an openstudio.Point3d."""
        return openstudio.Point3d(self.x, self.y, self.z)

    def __eq__(self, other):
        """Operator equal. Raises if not passed a Vertex."""
        if not isinstance(other, Vertex):
            raise NotImplementedError("Not implemented for any other types than Vertex itself")
        return isAlmostEqual3dPt(self, other)

    def __ne__(self, other):
        """Operator not equal."""
        return not self == other

    def __repr__(self):
        """Repr."""
        # return f"({self.x}, {self.y}, {self.z})"
        return f"({self.x:+.4f}, {self.y:+.4f}, {self.z:+.4f})"

__add__(other)

Return a + b.

Source code in geomeffibem/vertex.py 
 98
 99
100
def __add__(self, other) -> Vertex:
    """Return a + b."""
    return Vertex(x=self.x + other.x, y=self.y + other.y, z=self.z + other.z)

__eq__(other)

Operator equal. Raises if not passed a Vertex.

Source code in geomeffibem/vertex.py 
147
148
149
150
151
def __eq__(self, other):
    """Operator equal. Raises if not passed a Vertex."""
    if not isinstance(other, Vertex):
        raise NotImplementedError("Not implemented for any other types than Vertex itself")
    return isAlmostEqual3dPt(self, other)

__floordiv__(other)

Divides each coordinate by a scalar.

Source code in geomeffibem/vertex.py 
133
134
135
136
137
def __floordiv__(self, other) -> Vertex:
    """Divides each coordinate by a scalar."""
    if not isinstance(other, int) and not isinstance(other, float):
        raise ValueError("Division of a vertex by something else than a numeric is not supported")
    return Vertex(self.x // other, self.y // other, self.z // other)

__init__(x, y, z)

Vertex constructor.

Source code in geomeffibem/vertex.py 
32
33
34
35
36
37
def __init__(self, x, y, z):
    """Vertex constructor."""
    self.x = x
    self.y = y
    self.z = z
    self.surface = None

__mul__(other)

Multiplies each coordinate by a scalar.

Source code in geomeffibem/vertex.py 
115
116
117
118
119
def __mul__(self, other) -> Vertex:
    """Multiplies each coordinate by a scalar."""
    if not isinstance(other, int) and not isinstance(other, float):
        raise ValueError("Multiplication of a vertex by something else than a numeric is not supported")
    return Vertex(self.x * other, self.y * other, self.z * other)

__ne__(other)

Operator not equal.

Source code in geomeffibem/vertex.py 
153
154
155
def __ne__(self, other):
    """Operator not equal."""
    return not self == other

__neg__()

Return obj negated (-obj).

Source code in geomeffibem/vertex.py 
102
103
104
def __neg__(self):
    """Return obj negated (-obj)."""
    return Vertex(-self.x, -self.y, -self.z)

__repr__()

Repr.

Source code in geomeffibem/vertex.py 
157
158
159
160
def __repr__(self):
    """Repr."""
    # return f"({self.x}, {self.y}, {self.z})"
    return f"({self.x:+.4f}, {self.y:+.4f}, {self.z:+.4f})"

__rmul__(other)

Multiplies each coordinate by a scalar.

Source code in geomeffibem/vertex.py 
121
122
123
124
125
def __rmul__(self, other) -> Vertex:
    """Multiplies each coordinate by a scalar."""
    if not isinstance(other, int) and not isinstance(other, float):
        raise ValueError("Multiplication of a vertex by something else than a numeric is not supported")
    return self.__mul__(other)

__sub__(other)

Return a - b.

Source code in geomeffibem/vertex.py 
106
107
108
def __sub__(self, other) -> Vertex:
    """Return a - b."""
    return self + -other

__truediv__(other)

Divides each coordinate by a scalar.

Source code in geomeffibem/vertex.py 
127
128
129
130
131
def __truediv__(self, other) -> Vertex:
    """Divides each coordinate by a scalar."""
    if not isinstance(other, int) and not isinstance(other, float):
        raise ValueError("Division of a vertex by something else than a numeric is not supported")
    return Vertex(self.x / other, self.y / other, self.z / other)

copy()

Make a copy of this Vertex.

Source code in geomeffibem/vertex.py 
39
40
41
def copy(self):
    """Make a copy of this Vertex."""
    return Vertex(self.x, self.y, self.z)

cross(other, normalize=False)

Computes the cross product (a x b), which is a vector perpendicular to both a and b.

Source code in geomeffibem/vertex.py 
83
84
85
86
87
88
89
90
91
92
def cross(self, other, normalize: bool = False) -> Vertex:
    """Computes the cross product (a x b), which is a vector perpendicular to both a and b."""
    v = Vertex(
        x=(self.y * other.z - self.z * other.y),
        y=(self.z * other.x - self.x * other.z),
        z=(self.x * other.y - self.y * other.x),
    )
    if normalize:
        return v.normalize()
    return v

dot(other)

Computes the dot / scalar / inner product product of two vectors.

(a.b).

Source code in geomeffibem/vertex.py 
75
76
77
78
79
80
81
def dot(self, other) -> float:
    """Computes the dot / scalar / inner product product of two vectors.

    (a.b).
    """
    # return np.dot(v.to_numpy(), v2.to_numpy())
    return self.x * other.x + self.y * other.y + self.z * other.z

from_Point3d(pt) staticmethod

Factory method to construct from an openstudio.Point3d.

Source code in geomeffibem/vertex.py 
27
28
29
30
@staticmethod
def from_Point3d(pt: openstudio.Point3d):
    """Factory method to construct from an openstudio.Point3d."""
    return Vertex(pt.x(), pt.y(), pt.z())

from_numpy(arr) staticmethod

Factory method to construct from a numpy array (or a list) of 3 coords.

Source code in geomeffibem/vertex.py 
18
19
20
21
22
23
24
25
@staticmethod
def from_numpy(arr):
    """Factory method to construct from a numpy array (or a list) of 3 coords."""
    if not isinstance(arr, np.ndarray):
        arr = np.array(arr)
    if arr.shape != (3,):
        raise ValueError(f"Expected a numpy array with a dimension (3, ) (a Vector3d), got {arr.shape}")
    return Vertex(arr[0], arr[1], arr[2])

get_coords_on_plane(plane='xy')

Returns two coordinates on a given plane.

Source code in geomeffibem/vertex.py 
43
44
45
46
47
48
49
50
51
52
def get_coords_on_plane(self, plane='xy') -> Tuple[float, float]:
    """Returns two coordinates on a given plane."""
    if plane == 'xy':
        return self.x, self.y
    elif plane == 'xz':
        return self.x, self.z
    elif plane == 'yz':
        return self.y, self.z

    raise ValueError("Expected plane to be 'xy', 'xz' or 'yz'.")

length()

Get the length of the vector.

Source code in geomeffibem/vertex.py 
54
55
56
def length(self) -> float:
    """Get the length of the vector."""
    return np.sqrt(np.sum(self.to_numpy() ** 2))

normalize()

Normalize to a length of 1, returns a copy.

Source code in geomeffibem/vertex.py 
58
59
60
61
62
def normalize(self) -> Vertex:
    """Normalize to a length of 1, returns a copy."""
    v = self.copy()
    v.setLength(1.0)
    return v

outer_product(other)

Compute the outer product of this by another vector.

Source code in geomeffibem/vertex.py 
94
95
96
def outer_product(self, other) -> np.ndarray:
    """Compute the outer product of this by another vector."""
    return np.outer(self.to_numpy(), other.to_numpy())

setLength(newLength)

Change length of vector, in place.

Source code in geomeffibem/vertex.py 
64
65
66
67
68
69
70
71
72
73
def setLength(self, newLength: float) -> None:
    """Change length of vector, in place."""
    currentLength = self.length()
    if currentLength > 0:
        mult = newLength / currentLength
        self.x *= mult
        self.y *= mult
        self.z *= mult
    else:
        raise ValueError("Cannot normalize a vector of length 0")

to_Point3d()

Export to an openstudio.Point3d.

Source code in geomeffibem/vertex.py 
143
144
145
def to_Point3d(self) -> openstudio.Point3d:
    """Export to an openstudio.Point3d."""
    return openstudio.Point3d(self.x, self.y, self.z)

to_numpy()

Export to a numpy array of 3 coordinates.

Source code in geomeffibem/vertex.py 
139
140
141
def to_numpy(self) -> np.ndarray:
    """Export to a numpy array of 3 coordinates."""
    return np.array([self.x, self.y, self.z])

distance(lhs, rhs)

Distance between two vertices.

Source code in geomeffibem/vertex.py 
172
173
174
175
176
def distance(lhs: Vertex, rhs: Vertex) -> float:
    """Distance between two vertices."""
    squared_dist = np.sum((lhs.to_numpy() - rhs.to_numpy()) ** 2, axis=0)
    dist = np.sqrt(squared_dist)
    return dist

distanceFromPointToLine(start, end, test)

Distance between a point and a line.

Source code in geomeffibem/vertex.py 
179
180
181
182
183
184
def distanceFromPointToLine(start: Vertex, end: Vertex, test: Vertex) -> np.floating[Any]:
    """Distance between a point and a line."""
    s = start.to_numpy()
    e = end.to_numpy()
    p = test.to_numpy()
    return np.linalg.norm(np.cross(e - s, p - s) / np.linalg.norm(e - s))

getAngle(start, end)

Returns the angle between two vectors, in radians.

Source code in geomeffibem/vertex.py 
198
199
200
201
202
def getAngle(start: Vertex, end: Vertex) -> float:
    """Returns the angle between two vectors, in radians."""
    start = start.normalize()
    end = end.normalize()
    return np.arccos(start.dot(end))

getNewellVector(points)

Compute Newell vector from a list of points, direction is same as outward normal magnitude is twice the area.

Source code in geomeffibem/vertex.py 
205
206
207
208
209
210
211
212
213
214
215
216
217
def getNewellVector(points: List[Vertex]) -> Vertex:
    """Compute Newell vector from a list of points, direction is same as outward normal magnitude is twice the area."""
    n = len(points)
    if n < 3:
        raise ValueError("Cannot compute Newell Vector for less than 3 points")

    newellVector = Vertex(x=0, y=0, z=0)
    for i in range(n - 1):
        v1 = points[i] - points[0]
        v2 = points[i + 1] - points[0]
        newellVector += v1.cross(v2)

    return newellVector

getOutwardNormal(points)

Compute outward normal from a list of points.

Source code in geomeffibem/vertex.py 
220
221
222
223
def getOutwardNormal(points: list[Vertex]) -> Vertex:
    """Compute outward normal from a list of points."""
    newellVector = getNewellVector(points)
    return newellVector.normalize()

isAlmostEqual3dPt(v1, v2, tol=0.0127)

Checks if both vertices almost equal within tolerance.

Source code in geomeffibem/vertex.py 
166
167
168
169
def isAlmostEqual3dPt(v1: Vertex, v2: Vertex, tol=0.0127) -> bool:
    """Checks if both vertices almost equal within tolerance."""
    # 0.0127 m = 1.27 cm = 1/2 inch
    return not (abs((v1.to_numpy() - v2.to_numpy())) >= tol).any()

isPointOnLineBetweenPoints(start, end, test, tol=0.0127)

Checks whether a Vertex is on a segment.

If the distance(start, test) + distance(test, end) == distance(start, end) then it's on a line But we first check that the distance from the point to the line is also inferior to this tolerance

Source code in geomeffibem/vertex.py 
187
188
189
190
191
192
193
194
195
def isPointOnLineBetweenPoints(start: Vertex, end: Vertex, test: Vertex, tol: float = 0.0127) -> bool:
    """Checks whether a Vertex is on a segment.

    If the distance(start, test) + distance(test, end) == distance(start, end) then it's on a line
    But we first check that the distance from the point to the line is also inferior to this tolerance
    """
    if distanceFromPointToLine(start=start, end=end, test=test) < tol:
        return abs(distance(start, end) - (distance(start, test) + distance(test, end))) < tol
    return False

Surface

Surface and Surface3dEdge objects.

A Surface is a collection of Vertex. A Surface3dEdge is a side of a Surface.

Surface

A 3D Surface.

Source code in geomeffibem/surface.py 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
class Surface:
    """A 3D Surface."""

    @staticmethod
    def from_numpy_array(arr) -> Surface:
        """Factory method to construct from a numpy array of 3-coordinates arrays."""
        if isinstance(arr, list):
            arr = np.array(arr)
        if arr.shape[0] < 3:
            raise ValueError("Need at least 3 vertices to construct a Surface")
        if arr.shape[1] != 3:
            raise ValueError(f"Expected a numpy array with a dimension (N, 3), got {arr.shape}")
        return Surface([Vertex.from_numpy(x) for x in arr])

    @staticmethod
    def from_Point3dVector(points: Union[openstudio.Point3dVector, List[openstudio.Point3d]]) -> Surface:
        """Factory method to construct from an openstudio Point3dVector or a list of Point3d."""
        return Surface(vertices=[Vertex.from_Point3d(x) for x in points])

    @staticmethod
    def from_Surface(openstudio_surface: openstudio.model.Surface) -> Surface:
        """Factory method to construct from an openstudio.model.Surface."""
        if not isinstance(openstudio_surface, openstudio.model.Surface):
            raise ValueError("Expected an openstudio.model.Surface")
        return Surface(
            vertices=[Vertex.from_Point3d(x) for x in openstudio_surface.vertices()],
            name=openstudio_surface.nameString(),
        )

    @staticmethod
    def Floor(min_x=0.0, max_x=10.0, min_y=0.0, max_y=10.0, z=0.0) -> Surface:
        """Create a rectangular floor Surface (outward normal pointing down)."""
        # Counterclockwise, ULC convention, except here we want to create a floor so
        # outward normal must be pointing DOWN, so clockwise order
        vertices_arr = np.array(
            [
                [max_x, max_y, z],
                [max_x, min_y, z],
                [min_x, min_y, z],
                [min_x, max_y, z],
            ]
        )
        vertices = [Vertex.from_numpy(x) for x in vertices_arr]

        return Surface(vertices=vertices)

    @staticmethod
    def Rectangle(min_x=0.0, max_x=10.0, min_y=0.0, max_y=10.0, min_z=0.0, max_z=0.0) -> Surface:
        """Factory method to easily create a rectangular Surface, with ULC convention."""
        if abs(max_z - min_z) < 0.01:
            z = min_z
            if abs(z) < 0.01:
                print(
                    "Looks like you're trying to create a Floor surface... "
                    "use the Surface.Floor factory method if that's the case so outwardNormal points down."
                )

            vertices_arr = np.array(
                [
                    [min_x, max_y, z],  # Upper Left Corner
                    [min_x, min_y, z],  # Lower Left Corner
                    [max_x, min_y, z],  # Lower Right Corner
                    [max_x, max_y, z],  # Upper Right Corner
                ]
            )
        elif abs(max_x - min_x) < 0.01:
            x = min_x
            vertices_arr = np.array(
                [
                    [x, min_y, max_z],  # Upper Left Corner
                    [x, min_y, min_z],  # Lower Left Corner
                    [x, max_y, min_z],  # Lower Right Corner
                    [x, max_y, max_z],  # Upper Right Corner
                ]
            )

        elif abs(max_y - max_y) < 0.01:
            y = min_y
            vertices_arr = np.array(
                [
                    [min_x, y, max_z],  # Upper Left Corner
                    [min_x, y, min_z],  # Lower Left Corner
                    [max_x, y, min_z],  # Lower Right Corne
                    [max_x, y, max_z],  # Upper Right Corner
                ]
            )
        else:
            print("We expected at least one of x, y, z to be fixed, results are not guaranteed to work")
            vertices_arr = np.array(
                [
                    [min_x, min_y, max_z],
                    [min_x, max_y, min_z],
                    [max_x, max_y, min_z],
                    [max_x, min_y, max_z],
                ]
            )

        vertices = [Vertex.from_numpy(x) for x in vertices_arr]

        return Surface(vertices=vertices)

    def __init__(self, vertices, name=None):
        """Surface constructor."""
        if not isinstance(vertices, np.ndarray) and not isinstance(vertices, list):
            raise ValueError("Expected a list or numpy array of Vertex")

        for i, vertex in enumerate(vertices):
            if not isinstance(vertex, Vertex):
                raise ValueError(f"Element {i} is not a Vertex object")

        self.name = name
        self.os_plane = None

        self.vertices = copy.deepcopy(vertices)
        for vertex in self.vertices:
            vertex.surface = self

    def get_plane(self) -> Plane:
        """Returns the Plane of the Surface."""
        if self.os_plane is not None:
            return self.os_plane
        plane = openstudio.Plane(self.to_Point3dVector())
        self.os_plane = Plane(plane.a(), plane.b(), plane.c(), plane.d())
        return self.os_plane

    def get_plot_axis(self) -> str:
        """Returns a string representation of the plane it is on.

        TODO: raises if not exactly on 'xy', 'xz' or 'yz'
        """
        plane = self.get_plane()
        tol = 0.001
        if abs(abs(plane.a) - 1) < tol:
            return 'yz'
        if abs(abs(plane.b) - 1) < tol:
            return 'xz'
        if abs(abs(plane.c) - 1) < tol:
            return 'xy'

        # TODO
        raise NotImplementedError("Surface is not on a standard plane!")

    def plane(self) -> Plane:
        """Compute the plane from outwardNormal and the first point, not using OpenStudio."""
        normalVector = self.outwardNormal()
        if not np.isclose(normalVector.length(), 1.0):
            raise ValueError("Normal Unit Vector doesn't appear to be a unit vector")
        self.vertices[0]

        # d = -normalVector.x() * point.x() - normalVector.y() * point.y() - normalVector.z() * point.z();
        d = (-normalVector).dot(self.vertices[0])

        p = Plane(normalVector.x, normalVector.y, normalVector.z, d)
        for i, v in enumerate(self.vertices):
            if not p.pointOnPlane(v):
                print(f"Vertex {i} is not on the plane")
        return p

    def area(self) -> float:
        """Compute area of the surface."""
        newellVector = getNewellVector(self.vertices)
        return newellVector.length() / 2.0

    def outwardNormal(self) -> Vertex:
        """Returns the outward normal (normal unit vector)."""
        return getOutwardNormal(self.vertices)

    def tilt(self) -> float:
        """Returns the tilt of the surface, in radians, that is the angle between the outwardNormal and the Z axis."""
        z = Vertex(0.0, 0.0, 1.0)
        return getAngle(self.outwardNormal(), z)

    def azimuth(self) -> float:
        """Returns the azimuth of the surface, in radians.

        That is the angle between the outwardNormal and the North axis (Y-axis).
        """
        normal = self.outwardNormal()
        north = Vertex(0.0, 1.0, 0.0)
        angle = getAngle(normal, north)
        if normal.x < 0:
            return -angle + 2.0 * np.pi
        return angle

    def os_area(self) -> Vertex:
        """Returns area of the surface via openstudio."""
        return openstudio.getArea(self.to_Point3dVector()).get()

    def perimeter(self) -> float:
        """Returns the perimeter of the surface."""
        return sum([edge.length() for edge in self.to_Surface3dEdges()])

    def rough_centroid(self) -> Vertex:
        """Returns the centroid calculated in a rough way: the mean of the coordinates."""
        return Vertex.from_numpy(np.array([x.to_numpy() for x in self.vertices]).mean(axis=0))

    def os_centroid(self) -> Vertex:
        """Returns the centroid via openstudio."""
        centroid_ = openstudio.getCentroid(self.to_Point3dVector())
        if not centroid_.is_initialized():
            raise ValueError("OpenStudio failed to calculate centroid")
        return Vertex.from_Point3d(centroid_.get())

    def to_Point3dVector(self) -> List[openstudio.Point3d]:
        """Converts vertices to a list openstudio.Point3d."""
        return [v.to_Point3d() for v in self.vertices]

    def to_OSSurface(self, model: openstudio.model.Model) -> openstudio.model.Surface:
        """Creates an openstudio.model.Surface in the model passed as argument."""
        return openstudio.model.Surface(self.to_Point3dVector(), model)

    def to_numpy(self) -> np.ndarray:
        """Get a numpy array representing the vertices."""
        return np.array([v.to_numpy() for v in self.vertices])

    def to_Surface3dEdges(self) -> List[Surface3dEge]:
        """Converts vertex pairs to Surface3dEge."""
        edges = []
        for i, curVertex in enumerate(self.vertices):
            if i == len(self.vertices) - 1:
                nextVertex = self.vertices[0]
            else:
                nextVertex = self.vertices[i + 1]
            edges.append(Surface3dEge(start=curVertex, end=nextVertex, firstSurface=self))
        return edges

    def split_into_n_segments(self, n_segments, axis=None, plot=False) -> List[Surface]:
        """Splits a surface in N equal segments.

        If axis is not passed, it defaults to the first one of the plane
        eg: for a plane 'xy' it splits on 'x'
        """
        plot_axis = self.get_plot_axis()
        if axis is None:
            axis = plot_axis[0]
        if axis not in plot_axis:
            raise ValueError(f"This surface's plane is '{plot_axis}', so can't split on {axis=}")

        if n_segments < 2:
            raise ValueError("At least 2 segments needed")

        axis_to_index = {'x': 0, 'y': 1, 'z': 2}
        idx = axis_to_index[axis]
        v_np = self.to_numpy()
        minimum = v_np[:, idx].min()
        maximum = v_np[:, idx].max()
        segment_length = (maximum - minimum) / n_segments
        is_max = v_np[:, idx] == maximum
        is_min = ~is_max

        cur_min = minimum
        cur_max = cur_min + segment_length

        new_surfaces = []

        for i in range(n_segments):
            # print(cur_min, cur_max)
            v_np_i = v_np.copy()
            v_np_i[is_min, idx] = cur_min
            v_np_i[is_max, idx] = cur_max
            new_surface = Surface.from_numpy_array(v_np_i)
            if self.name:
                new_surface.name = f'{self.name}-{i+1}'
            new_surfaces.append(new_surface)

            cur_min = cur_max
            cur_max += segment_length

        if plot:
            fig, ax = plt.subplots(figsize=(16, 9))
            for new_surface in new_surfaces:
                new_surface.plot(ax=ax)

        return new_surfaces

    def rotate(self, degrees: float, axis=None) -> Surface:
        """Rotates a surface by an amount of degrees.

        Args:
        -----
        * degrees (float): the angle to rotate it by, in degrees. Positive means clockwise
        * axis (Vertex): if none, uses the Z axis

        Returns:
        ---------
        * a new Surface object with rotated vertices
        """
        if axis is None:
            axis = Vertex(0.0, 0.0, 1.0)

        # Lazy load to avoid circular import
        from geomeffibem.transformation import Transformation

        return Transformation.Rotation(axis=axis, radians=-openstudio.degToRad(degrees)) * self

    def translate(self, translation: Vertex) -> Surface:
        """Translates a surface along a translation vector."""
        from geomeffibem.transformation import Transformation

        return Transformation.Translation(translation=translation) * self

    def plot(self, name: Union[bool, str] = True, **kwargs):
        """Calls plot_vertices, cf help(plot_vertices)."""
        if isinstance(name, str):
            name = name
        elif name:
            name = self.name
        return plot_vertices(surface_like=self, name=name, **kwargs)

    def __repr__(self):
        """Repr."""
        s = ""
        if self.name is not None:
            s += f"Surface '{self.name}' = "
        s += "["
        if self.name is not None:
            s += "\n "
        imax = len(self.vertices) - 1
        for i, v in enumerate(self.vertices):
            if i > 0:
                s += " "
            s += f"{v}"
            if i < imax:
                s += ",\n"
        s += "]"
        return s

Floor(min_x=0.0, max_x=10.0, min_y=0.0, max_y=10.0, z=0.0) staticmethod

Create a rectangular floor Surface (outward normal pointing down).

Source code in geomeffibem/surface.py 
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
@staticmethod
def Floor(min_x=0.0, max_x=10.0, min_y=0.0, max_y=10.0, z=0.0) -> Surface:
    """Create a rectangular floor Surface (outward normal pointing down)."""
    # Counterclockwise, ULC convention, except here we want to create a floor so
    # outward normal must be pointing DOWN, so clockwise order
    vertices_arr = np.array(
        [
            [max_x, max_y, z],
            [max_x, min_y, z],
            [min_x, min_y, z],
            [min_x, max_y, z],
        ]
    )
    vertices = [Vertex.from_numpy(x) for x in vertices_arr]

    return Surface(vertices=vertices)

Rectangle(min_x=0.0, max_x=10.0, min_y=0.0, max_y=10.0, min_z=0.0, max_z=0.0) staticmethod

Factory method to easily create a rectangular Surface, with ULC convention.

Source code in geomeffibem/surface.py 
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
@staticmethod
def Rectangle(min_x=0.0, max_x=10.0, min_y=0.0, max_y=10.0, min_z=0.0, max_z=0.0) -> Surface:
    """Factory method to easily create a rectangular Surface, with ULC convention."""
    if abs(max_z - min_z) < 0.01:
        z = min_z
        if abs(z) < 0.01:
            print(
                "Looks like you're trying to create a Floor surface... "
                "use the Surface.Floor factory method if that's the case so outwardNormal points down."
            )

        vertices_arr = np.array(
            [
                [min_x, max_y, z],  # Upper Left Corner
                [min_x, min_y, z],  # Lower Left Corner
                [max_x, min_y, z],  # Lower Right Corner
                [max_x, max_y, z],  # Upper Right Corner
            ]
        )
    elif abs(max_x - min_x) < 0.01:
        x = min_x
        vertices_arr = np.array(
            [
                [x, min_y, max_z],  # Upper Left Corner
                [x, min_y, min_z],  # Lower Left Corner
                [x, max_y, min_z],  # Lower Right Corner
                [x, max_y, max_z],  # Upper Right Corner
            ]
        )

    elif abs(max_y - max_y) < 0.01:
        y = min_y
        vertices_arr = np.array(
            [
                [min_x, y, max_z],  # Upper Left Corner
                [min_x, y, min_z],  # Lower Left Corner
                [max_x, y, min_z],  # Lower Right Corne
                [max_x, y, max_z],  # Upper Right Corner
            ]
        )
    else:
        print("We expected at least one of x, y, z to be fixed, results are not guaranteed to work")
        vertices_arr = np.array(
            [
                [min_x, min_y, max_z],
                [min_x, max_y, min_z],
                [max_x, max_y, min_z],
                [max_x, min_y, max_z],
            ]
        )

    vertices = [Vertex.from_numpy(x) for x in vertices_arr]

    return Surface(vertices=vertices)

__init__(vertices, name=None)

Surface constructor.

Source code in geomeffibem/surface.py 
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
def __init__(self, vertices, name=None):
    """Surface constructor."""
    if not isinstance(vertices, np.ndarray) and not isinstance(vertices, list):
        raise ValueError("Expected a list or numpy array of Vertex")

    for i, vertex in enumerate(vertices):
        if not isinstance(vertex, Vertex):
            raise ValueError(f"Element {i} is not a Vertex object")

    self.name = name
    self.os_plane = None

    self.vertices = copy.deepcopy(vertices)
    for vertex in self.vertices:
        vertex.surface = self

__repr__()

Repr.

Source code in geomeffibem/surface.py 
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
def __repr__(self):
    """Repr."""
    s = ""
    if self.name is not None:
        s += f"Surface '{self.name}' = "
    s += "["
    if self.name is not None:
        s += "\n "
    imax = len(self.vertices) - 1
    for i, v in enumerate(self.vertices):
        if i > 0:
            s += " "
        s += f"{v}"
        if i < imax:
            s += ",\n"
    s += "]"
    return s

area()

Compute area of the surface.

Source code in geomeffibem/surface.py 
242
243
244
245
def area(self) -> float:
    """Compute area of the surface."""
    newellVector = getNewellVector(self.vertices)
    return newellVector.length() / 2.0

azimuth()

Returns the azimuth of the surface, in radians.

That is the angle between the outwardNormal and the North axis (Y-axis).

Source code in geomeffibem/surface.py 
256
257
258
259
260
261
262
263
264
265
266
def azimuth(self) -> float:
    """Returns the azimuth of the surface, in radians.

    That is the angle between the outwardNormal and the North axis (Y-axis).
    """
    normal = self.outwardNormal()
    north = Vertex(0.0, 1.0, 0.0)
    angle = getAngle(normal, north)
    if normal.x < 0:
        return -angle + 2.0 * np.pi
    return angle

from_Point3dVector(points) staticmethod

Factory method to construct from an openstudio Point3dVector or a list of Point3d.

Source code in geomeffibem/surface.py 
 98
 99
100
101
@staticmethod
def from_Point3dVector(points: Union[openstudio.Point3dVector, List[openstudio.Point3d]]) -> Surface:
    """Factory method to construct from an openstudio Point3dVector or a list of Point3d."""
    return Surface(vertices=[Vertex.from_Point3d(x) for x in points])

from_Surface(openstudio_surface) staticmethod

Factory method to construct from an openstudio.model.Surface.

Source code in geomeffibem/surface.py 
103
104
105
106
107
108
109
110
111
@staticmethod
def from_Surface(openstudio_surface: openstudio.model.Surface) -> Surface:
    """Factory method to construct from an openstudio.model.Surface."""
    if not isinstance(openstudio_surface, openstudio.model.Surface):
        raise ValueError("Expected an openstudio.model.Surface")
    return Surface(
        vertices=[Vertex.from_Point3d(x) for x in openstudio_surface.vertices()],
        name=openstudio_surface.nameString(),
    )

from_numpy_array(arr) staticmethod

Factory method to construct from a numpy array of 3-coordinates arrays.

Source code in geomeffibem/surface.py 
87
88
89
90
91
92
93
94
95
96
@staticmethod
def from_numpy_array(arr) -> Surface:
    """Factory method to construct from a numpy array of 3-coordinates arrays."""
    if isinstance(arr, list):
        arr = np.array(arr)
    if arr.shape[0] < 3:
        raise ValueError("Need at least 3 vertices to construct a Surface")
    if arr.shape[1] != 3:
        raise ValueError(f"Expected a numpy array with a dimension (N, 3), got {arr.shape}")
    return Surface([Vertex.from_numpy(x) for x in arr])

get_plane()

Returns the Plane of the Surface.

Source code in geomeffibem/surface.py 
201
202
203
204
205
206
207
def get_plane(self) -> Plane:
    """Returns the Plane of the Surface."""
    if self.os_plane is not None:
        return self.os_plane
    plane = openstudio.Plane(self.to_Point3dVector())
    self.os_plane = Plane(plane.a(), plane.b(), plane.c(), plane.d())
    return self.os_plane

get_plot_axis()

Returns a string representation of the plane it is on.

TODO: raises if not exactly on 'xy', 'xz' or 'yz'

Source code in geomeffibem/surface.py 
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
def get_plot_axis(self) -> str:
    """Returns a string representation of the plane it is on.

    TODO: raises if not exactly on 'xy', 'xz' or 'yz'
    """
    plane = self.get_plane()
    tol = 0.001
    if abs(abs(plane.a) - 1) < tol:
        return 'yz'
    if abs(abs(plane.b) - 1) < tol:
        return 'xz'
    if abs(abs(plane.c) - 1) < tol:
        return 'xy'

    # TODO
    raise NotImplementedError("Surface is not on a standard plane!")

os_area()

Returns area of the surface via openstudio.

Source code in geomeffibem/surface.py 
268
269
270
def os_area(self) -> Vertex:
    """Returns area of the surface via openstudio."""
    return openstudio.getArea(self.to_Point3dVector()).get()

os_centroid()

Returns the centroid via openstudio.

Source code in geomeffibem/surface.py 
280
281
282
283
284
285
def os_centroid(self) -> Vertex:
    """Returns the centroid via openstudio."""
    centroid_ = openstudio.getCentroid(self.to_Point3dVector())
    if not centroid_.is_initialized():
        raise ValueError("OpenStudio failed to calculate centroid")
    return Vertex.from_Point3d(centroid_.get())

outwardNormal()

Returns the outward normal (normal unit vector).

Source code in geomeffibem/surface.py 
247
248
249
def outwardNormal(self) -> Vertex:
    """Returns the outward normal (normal unit vector)."""
    return getOutwardNormal(self.vertices)

perimeter()

Returns the perimeter of the surface.

Source code in geomeffibem/surface.py 
272
273
274
def perimeter(self) -> float:
    """Returns the perimeter of the surface."""
    return sum([edge.length() for edge in self.to_Surface3dEdges()])

plane()

Compute the plane from outwardNormal and the first point, not using OpenStudio.

Source code in geomeffibem/surface.py 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
def plane(self) -> Plane:
    """Compute the plane from outwardNormal and the first point, not using OpenStudio."""
    normalVector = self.outwardNormal()
    if not np.isclose(normalVector.length(), 1.0):
        raise ValueError("Normal Unit Vector doesn't appear to be a unit vector")
    self.vertices[0]

    # d = -normalVector.x() * point.x() - normalVector.y() * point.y() - normalVector.z() * point.z();
    d = (-normalVector).dot(self.vertices[0])

    p = Plane(normalVector.x, normalVector.y, normalVector.z, d)
    for i, v in enumerate(self.vertices):
        if not p.pointOnPlane(v):
            print(f"Vertex {i} is not on the plane")
    return p

plot(name=True, **kwargs)

Calls plot_vertices, cf help(plot_vertices).

Source code in geomeffibem/surface.py 
385
386
387
388
389
390
391
def plot(self, name: Union[bool, str] = True, **kwargs):
    """Calls plot_vertices, cf help(plot_vertices)."""
    if isinstance(name, str):
        name = name
    elif name:
        name = self.name
    return plot_vertices(surface_like=self, name=name, **kwargs)

rotate(degrees, axis=None)

Rotates a surface by an amount of degrees.

Args:
  • degrees (float): the angle to rotate it by, in degrees. Positive means clockwise
  • axis (Vertex): if none, uses the Z axis
Returns:
  • a new Surface object with rotated vertices
Source code in geomeffibem/surface.py 
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
def rotate(self, degrees: float, axis=None) -> Surface:
    """Rotates a surface by an amount of degrees.

    Args:
    -----
    * degrees (float): the angle to rotate it by, in degrees. Positive means clockwise
    * axis (Vertex): if none, uses the Z axis

    Returns:
    ---------
    * a new Surface object with rotated vertices
    """
    if axis is None:
        axis = Vertex(0.0, 0.0, 1.0)

    # Lazy load to avoid circular import
    from geomeffibem.transformation import Transformation

    return Transformation.Rotation(axis=axis, radians=-openstudio.degToRad(degrees)) * self

rough_centroid()

Returns the centroid calculated in a rough way: the mean of the coordinates.

Source code in geomeffibem/surface.py 
276
277
278
def rough_centroid(self) -> Vertex:
    """Returns the centroid calculated in a rough way: the mean of the coordinates."""
    return Vertex.from_numpy(np.array([x.to_numpy() for x in self.vertices]).mean(axis=0))

split_into_n_segments(n_segments, axis=None, plot=False)

Splits a surface in N equal segments.

If axis is not passed, it defaults to the first one of the plane eg: for a plane 'xy' it splits on 'x'

Source code in geomeffibem/surface.py 
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
def split_into_n_segments(self, n_segments, axis=None, plot=False) -> List[Surface]:
    """Splits a surface in N equal segments.

    If axis is not passed, it defaults to the first one of the plane
    eg: for a plane 'xy' it splits on 'x'
    """
    plot_axis = self.get_plot_axis()
    if axis is None:
        axis = plot_axis[0]
    if axis not in plot_axis:
        raise ValueError(f"This surface's plane is '{plot_axis}', so can't split on {axis=}")

    if n_segments < 2:
        raise ValueError("At least 2 segments needed")

    axis_to_index = {'x': 0, 'y': 1, 'z': 2}
    idx = axis_to_index[axis]
    v_np = self.to_numpy()
    minimum = v_np[:, idx].min()
    maximum = v_np[:, idx].max()
    segment_length = (maximum - minimum) / n_segments
    is_max = v_np[:, idx] == maximum
    is_min = ~is_max

    cur_min = minimum
    cur_max = cur_min + segment_length

    new_surfaces = []

    for i in range(n_segments):
        # print(cur_min, cur_max)
        v_np_i = v_np.copy()
        v_np_i[is_min, idx] = cur_min
        v_np_i[is_max, idx] = cur_max
        new_surface = Surface.from_numpy_array(v_np_i)
        if self.name:
            new_surface.name = f'{self.name}-{i+1}'
        new_surfaces.append(new_surface)

        cur_min = cur_max
        cur_max += segment_length

    if plot:
        fig, ax = plt.subplots(figsize=(16, 9))
        for new_surface in new_surfaces:
            new_surface.plot(ax=ax)

    return new_surfaces

tilt()

Returns the tilt of the surface, in radians, that is the angle between the outwardNormal and the Z axis.

Source code in geomeffibem/surface.py 
251
252
253
254
def tilt(self) -> float:
    """Returns the tilt of the surface, in radians, that is the angle between the outwardNormal and the Z axis."""
    z = Vertex(0.0, 0.0, 1.0)
    return getAngle(self.outwardNormal(), z)

to_OSSurface(model)

Creates an openstudio.model.Surface in the model passed as argument.

Source code in geomeffibem/surface.py 
291
292
293
def to_OSSurface(self, model: openstudio.model.Model) -> openstudio.model.Surface:
    """Creates an openstudio.model.Surface in the model passed as argument."""
    return openstudio.model.Surface(self.to_Point3dVector(), model)

to_Point3dVector()

Converts vertices to a list openstudio.Point3d.

Source code in geomeffibem/surface.py 
287
288
289
def to_Point3dVector(self) -> List[openstudio.Point3d]:
    """Converts vertices to a list openstudio.Point3d."""
    return [v.to_Point3d() for v in self.vertices]

to_Surface3dEdges()

Converts vertex pairs to Surface3dEge.

Source code in geomeffibem/surface.py 
299
300
301
302
303
304
305
306
307
308
def to_Surface3dEdges(self) -> List[Surface3dEge]:
    """Converts vertex pairs to Surface3dEge."""
    edges = []
    for i, curVertex in enumerate(self.vertices):
        if i == len(self.vertices) - 1:
            nextVertex = self.vertices[0]
        else:
            nextVertex = self.vertices[i + 1]
        edges.append(Surface3dEge(start=curVertex, end=nextVertex, firstSurface=self))
    return edges

to_numpy()

Get a numpy array representing the vertices.

Source code in geomeffibem/surface.py 
295
296
297
def to_numpy(self) -> np.ndarray:
    """Get a numpy array representing the vertices."""
    return np.array([v.to_numpy() for v in self.vertices])

translate(translation)

Translates a surface along a translation vector.

Source code in geomeffibem/surface.py 
379
380
381
382
383
def translate(self, translation: Vertex) -> Surface:
    """Translates a surface along a translation vector."""
    from geomeffibem.transformation import Transformation

    return Transformation.Translation(translation=translation) * self

Surface3dEge

An Edge has a start and an end Vertex, and a list of surfaces it was found on.

Source code in geomeffibem/surface.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
class Surface3dEge:
    """An Edge has a start and an end Vertex, and a list of surfaces it was found on."""

    def __init__(self, start: Vertex, end: Vertex, firstSurface: Surface):
        """Constructor."""
        self.start = start
        self.end = end
        self.allSurfaces = [firstSurface]

    def containsPoints(self, testVertex: Vertex) -> bool:
        """Checks whether a Point is on the edge.

        It is not almost equal to the start and end points, and,
        isPointOnLineBetweenPoints(start, end, testVertex) is true.
        """
        return (
            not isAlmostEqual3dPt(self.start, testVertex)
            and not isAlmostEqual3dPt(self.end, testVertex)
            and isPointOnLineBetweenPoints(self.start, self.end, testVertex)
        )

    def length(self) -> float:
        """Compute distance from start to end."""
        return distance(self.start, self.end)

    def count(self) -> int:
        """Number of Surfaces it was found on."""
        return len(self.allSurfaces)

    def __eq__(self, other):
        """Operator equal."""
        if not isinstance(other, Surface3dEge):
            raise NotImplementedError("Not implemented for any other types than Surface3dEge itself")

        return (isAlmostEqual3dPt(self.start, other.start) and isAlmostEqual3dPt(self.end, other.end)) or (
            isAlmostEqual3dPt(self.start, other.end) and isAlmostEqual3dPt(self.end, other.start)
        )

    def __ne__(self, other):
        """Operator not equal."""
        return not self == other

    def __repr__(self):
        """Repr."""
        return f"start={self.start}, end={self.end}, count={self.count()}, firstSurface={self.allSurfaces[0].name}"

    def plot_on_first_surface(self, ax=None):
        """Plots this segment in red atop the outline of the Surface it came from."""
        surface = self.allSurfaces[0]

        if ax is None:
            fig, ax = plt.subplots(figsize=(16, 9))
        surface.plot(ax=ax)
        plot_vertices([self.start, self.end], plane=surface.get_plot_axis(), c='r', ax=ax, annotate=False)

__eq__(other)

Operator equal.

Source code in geomeffibem/surface.py 
57
58
59
60
61
62
63
64
def __eq__(self, other):
    """Operator equal."""
    if not isinstance(other, Surface3dEge):
        raise NotImplementedError("Not implemented for any other types than Surface3dEge itself")

    return (isAlmostEqual3dPt(self.start, other.start) and isAlmostEqual3dPt(self.end, other.end)) or (
        isAlmostEqual3dPt(self.start, other.end) and isAlmostEqual3dPt(self.end, other.start)
    )

__init__(start, end, firstSurface)

Constructor.

Source code in geomeffibem/surface.py 
31
32
33
34
35
def __init__(self, start: Vertex, end: Vertex, firstSurface: Surface):
    """Constructor."""
    self.start = start
    self.end = end
    self.allSurfaces = [firstSurface]

__ne__(other)

Operator not equal.

Source code in geomeffibem/surface.py 
66
67
68
def __ne__(self, other):
    """Operator not equal."""
    return not self == other

__repr__()

Repr.

Source code in geomeffibem/surface.py 
70
71
72
def __repr__(self):
    """Repr."""
    return f"start={self.start}, end={self.end}, count={self.count()}, firstSurface={self.allSurfaces[0].name}"

containsPoints(testVertex)

Checks whether a Point is on the edge.

It is not almost equal to the start and end points, and, isPointOnLineBetweenPoints(start, end, testVertex) is true.

Source code in geomeffibem/surface.py 
37
38
39
40
41
42
43
44
45
46
47
def containsPoints(self, testVertex: Vertex) -> bool:
    """Checks whether a Point is on the edge.

    It is not almost equal to the start and end points, and,
    isPointOnLineBetweenPoints(start, end, testVertex) is true.
    """
    return (
        not isAlmostEqual3dPt(self.start, testVertex)
        and not isAlmostEqual3dPt(self.end, testVertex)
        and isPointOnLineBetweenPoints(self.start, self.end, testVertex)
    )

count()

Number of Surfaces it was found on.

Source code in geomeffibem/surface.py 
53
54
55
def count(self) -> int:
    """Number of Surfaces it was found on."""
    return len(self.allSurfaces)

length()

Compute distance from start to end.

Source code in geomeffibem/surface.py 
49
50
51
def length(self) -> float:
    """Compute distance from start to end."""
    return distance(self.start, self.end)

plot_on_first_surface(ax=None)

Plots this segment in red atop the outline of the Surface it came from.

Source code in geomeffibem/surface.py 
74
75
76
77
78
79
80
81
def plot_on_first_surface(self, ax=None):
    """Plots this segment in red atop the outline of the Surface it came from."""
    surface = self.allSurfaces[0]

    if ax is None:
        fig, ax = plt.subplots(figsize=(16, 9))
    surface.plot(ax=ax)
    plot_vertices([self.start, self.end], plane=surface.get_plot_axis(), c='r', ax=ax, annotate=False)

get_surface_from_surface_like(surface_like)

Helper to get a Surface (class) from a surface like object.

Source code in geomeffibem/surface.py 
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
def get_surface_from_surface_like(surface_like: Union[Surface, List[Vertex], openstudio.model.Surface]) -> Surface:
    """Helper to get a Surface (class) from a surface like object."""
    if isinstance(surface_like, openstudio.model.Surface):
        surface = Surface.from_Surface(surface_like)
    elif isinstance(surface_like, Surface):
        surface = surface_like
    else:
        if isinstance(surface_like, list):
            surface_like = np.array(surface_like)

        if isinstance(surface_like[0], openstudio.Point3d):
            surface = Surface.from_Point3dVector(surface_like)
        elif isinstance(surface_like[0], np.ndarray):
            surface = Surface.from_numpy_array(surface_like)
        elif isinstance(surface_like[0], Vertex):
            surface = Surface(surface_like)

    return surface

plot_vertices(surface_like, ax=None, center_axes=False, with_rough_centroid=False, with_os_centroid=False, annotate=True, linewidth=None, force_align=False, name=None, plane=None, annotate_kwargs=dict(color='r', xytext=(5, 5), textcoords='offset points'), **kwargs)

Plot any surface-like object in 2D.

Accepts a Surface, a list or numpy array of Vertex, or an openstudio.model.Surface object

TODO: Assumes the surface is planar and falls exactly on 'xy', 'xz' or 'yz' plane currently

Source code in geomeffibem/surface.py 
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
def plot_vertices(
    surface_like: Union[Surface, List[Vertex], openstudio.model.Surface],
    ax=None,
    center_axes=False,
    with_rough_centroid=False,
    with_os_centroid=False,
    annotate=True,
    linewidth=None,
    force_align=False,
    name=None,
    plane=None,
    annotate_kwargs=dict(color='r', xytext=(5, 5), textcoords='offset points'),
    # Passed to ax.plot/plt.plot
    **kwargs,
):
    """Plot any surface-like object in 2D.

    Accepts a Surface, a list or numpy array of Vertex, or an openstudio.model.Surface object

    TODO: Assumes the surface is planar and falls exactly on 'xy', 'xz' or 'yz' plane currently
    """
    surface = get_surface_from_surface_like(surface_like=surface_like)

    is_aligned = False
    if plane is None and not force_align:
        try:
            plane = surface.get_plot_axis()
        except NotImplementedError as e:
            print(e)
            force_align = True
    if plane is None and force_align:
        print("aligning Face")
        # Lazy load to avoid circular import
        from geomeffibem.transformation import Transformation

        faceTransformation = Transformation.alignFace(surface.vertices)
        faceTransformationInverse = faceTransformation.inverse()
        points = faceTransformationInverse * surface.vertices
        surface = Surface(points, name=surface.name)
        plane = 'xy'
        is_aligned = True

    points = surface.to_numpy()

    if plane == 'xy':
        xs = points[:, 0]
        ys = points[:, 1]
    elif plane == 'xz':
        xs = points[:, 0]
        ys = points[:, 2]
    elif plane == 'yz':
        xs = points[:, 1]
        ys = points[:, 2]
    else:
        raise ValueError("plane must be in ['xy', 'xz', 'yz']")
    if ax is None:
        # print("Making a figure")
        max_width = xs.max() - xs.min()
        max_height = ys.max() - ys.min()
        h = 6
        w = h * max_width / max_height
        # print(w, h)
        fig, ax = plt.subplots(figsize=(w, h))

    ax.plot(np.append(xs, xs[0]), np.append(ys, ys[0]), marker='x', markeredgecolor='r', linewidth=linewidth, **kwargs)
    ax.set_xlabel(plane[0] if not is_aligned else "x'")
    ax.set_ylabel(plane[1] if not is_aligned else "y'")
    if annotate:
        for i, (x, y) in enumerate(zip(xs, ys)):
            ax.annotate(f"{i+1} ({x:.2f}, {y:.2f})", xy=(x, y), **annotate_kwargs)

    if with_rough_centroid:
        centroid_x, centroid_y = surface.rough_centroid().get_coords_on_plane(plane=plane)
        ax.annotate(f"rough ({centroid_x}, {centroid_y})", xy=(centroid_x, centroid_y))
        ax.plot(centroid_x, centroid_y, 'rx')
    if with_os_centroid or name is not None and name is not False:
        centroid_x, centroid_y = surface.os_centroid().get_coords_on_plane(plane=plane)
        if with_os_centroid:
            ax.annotate(f"os ({centroid_x}, {centroid_y})", xy=(centroid_x, centroid_y))
            ax.plot(centroid_x, centroid_y, 'gx')
            if name:
                ax.annotate(
                    name + (" (aligned)" if is_aligned else ""),
                    xy=(centroid_x, centroid_y),
                    xytext=(0, 50),
                    textcoords='offset pixels',
                    color='b',
                    arrowprops=dict(edgecolor='b', lw=1, ls='-', arrowstyle='->'),
                )
        else:
            ax.annotate(
                name + (" (aligned)" if is_aligned else ""), xy=(centroid_x, centroid_y), ha='center', va='center'
            )

    ax.spines['right'].set_color('none')
    ax.spines['top'].set_color('none')  # hide top axis

    if center_axes:
        ax.spines['bottom'].set_position('zero')  # x-axis where y=0
        ax.spines['left'].set_position('zero')
        ax.xaxis.set_label_position("top")

    return ax

Plane

Plane.

Plane

A 3D Plane.

Equation is ax + by + cz + d = 0

Source code in geomeffibem/plane.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
class Plane:
    """A 3D Plane.

    Equation is `ax + by + cz + d = 0`
    """

    def __init__(self, a, b, c, d):
        """Plane constructor."""
        self.a = a
        self.b = b
        self.c = c
        self.d = d

    def outwardNormal(self) -> np.ndarray:
        """The outwardNormal of the plane."""
        return Vertex(self.a, self.b, self.c)

    def is_orthogonal(self) -> bool:
        """Checks if the plane is orthogonal."""
        return np.sum(np.abs(self.outwardNormal().to_numpy())) == 1.0

    def pointOnPlane(self, point: Vertex, tol=0.001) -> bool:
        """Checks whether the Vertex is on the Plane."""
        # project point to plane
        projected = self.project(point)

        return distance(point, projected) <= tol

    def project(self, point: Vertex) -> Vertex:
        """Project a point onto a Plane."""
        # http://www.9math.com/book/projection-point-plane
        if not isinstance(point, Vertex):
            raise ValueError("Expected a Vertex object")
        u = point.x
        v = point.y
        w = point.z

        num = self.a * u + self.b * v + self.c * w + self.d
        den = self.a * self.a + self.b * self.b + self.c * self.c  # this should always be 1.0
        ratio = num / den

        x = u - self.a * ratio
        y = v - self.b * ratio
        z = w - self.c * ratio

        return Vertex(x, y, z)

    def __repr__(self):
        """Repr."""
        return f"Plane ({self.a}, {self.b}, {self.c}, {self.d})"

__init__(a, b, c, d)

Plane constructor.

Source code in geomeffibem/plane.py 
14
15
16
17
18
19
def __init__(self, a, b, c, d):
    """Plane constructor."""
    self.a = a
    self.b = b
    self.c = c
    self.d = d

__repr__()

Repr.

Source code in geomeffibem/plane.py 
55
56
57
def __repr__(self):
    """Repr."""
    return f"Plane ({self.a}, {self.b}, {self.c}, {self.d})"

is_orthogonal()

Checks if the plane is orthogonal.

Source code in geomeffibem/plane.py 
25
26
27
def is_orthogonal(self) -> bool:
    """Checks if the plane is orthogonal."""
    return np.sum(np.abs(self.outwardNormal().to_numpy())) == 1.0

outwardNormal()

The outwardNormal of the plane.

Source code in geomeffibem/plane.py 
21
22
23
def outwardNormal(self) -> np.ndarray:
    """The outwardNormal of the plane."""
    return Vertex(self.a, self.b, self.c)

pointOnPlane(point, tol=0.001)

Checks whether the Vertex is on the Plane.

Source code in geomeffibem/plane.py 
29
30
31
32
33
34
def pointOnPlane(self, point: Vertex, tol=0.001) -> bool:
    """Checks whether the Vertex is on the Plane."""
    # project point to plane
    projected = self.project(point)

    return distance(point, projected) <= tol

project(point)

Project a point onto a Plane.

Source code in geomeffibem/plane.py 
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
def project(self, point: Vertex) -> Vertex:
    """Project a point onto a Plane."""
    # http://www.9math.com/book/projection-point-plane
    if not isinstance(point, Vertex):
        raise ValueError("Expected a Vertex object")
    u = point.x
    v = point.y
    w = point.z

    num = self.a * u + self.b * v + self.c * w + self.d
    den = self.a * self.a + self.b * self.b + self.c * self.c  # this should always be 1.0
    ratio = num / den

    x = u - self.a * ratio
    y = v - self.b * ratio
    z = w - self.c * ratio

    return Vertex(x, y, z)

Polyhedron

A Polyhedron is a collection of Surface objects.

It's meant to represent a volume. You can check whether it's enclosed or not.

Polyhedron

A collection of Surfaces, meant to represent a Volume.

Source code in geomeffibem/polyhedron.py 
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
class Polyhedron:
    """A collection of Surfaces, meant to represent a Volume."""

    def __init__(self, surfaces: List[Surface]):
        """Constructor from a list of Surface objects."""
        if not isinstance(surfaces, np.ndarray) and not isinstance(surfaces, list):
            raise ValueError("Expected a list or numpy array of Surfaces")

        for i, surface in enumerate(surfaces):
            if not isinstance(surface, Surface):
                raise ValueError(f"Element {i} is not a Surface object")
        self.surfaces = surfaces

    def get_surface_by_name(self, name):
        """Locate a surface by its name."""
        for s in self.surfaces:
            if s.name is not None and s.name == name:
                return s

    def numVertices(self):
        """Counts the total number of vertices for all surfaces."""
        count = 0
        for s in self.surfaces:
            count += len(s.vertices)
        return count

    def uniqueVertices(self) -> List[Vertex]:
        """Get a list of unique vertices (uses Vertex __eq__ operator which has a tolerance)."""
        uniqueVertices: List[Vertex] = []
        for s in self.surfaces:
            for vertex in s.vertices:
                found = False
                for unique_v in uniqueVertices:
                    if unique_v == vertex:
                        found = True
                        break
                if not found:
                    uniqueVertices.append(vertex)
        return uniqueVertices

    @staticmethod
    def edgesNotTwoForEnclosedVolumeTest(zonePoly: Polyhedron) -> Tuple[List[Surface3dEge], List[Surface3dEge]]:
        """Counts the number of times an Edge is used.

        Returns the ones that isn't used twice (and the ones used twice for debugging/inspection)
        """
        uniqueSurface3dEdges: List[Surface3dEge] = []

        for surface in zonePoly.surfaces:
            for edge in surface.to_Surface3dEdges():
                found = False
                for uniqEdge in uniqueSurface3dEdges:
                    if uniqEdge == edge:
                        uniqEdge.allSurfaces.append(surface)
                        found = True
                        break
                if not found:
                    uniqueSurface3dEdges.append(edge)

        edgesNotTwoCount = [x for x in uniqueSurface3dEdges if x.count() != 2]
        edgesTwoCount = [x for x in uniqueSurface3dEdges if x.count() == 2]
        return edgesNotTwoCount, edgesTwoCount

    def updateZonePolygonsForMissingColinearPoints(self) -> Polyhedron:
        """Creates a new Polyhedron with extra vertices when a point is found to be on a line segment."""
        updZonePoly = copy.deepcopy(self)

        uniqVertices = self.uniqueVertices()

        for surface in updZonePoly.surfaces:
            insertedVertex = True
            while insertedVertex:
                insertedVertex = False
                for i, edge in enumerate(surface.to_Surface3dEdges()):
                    for testVertex in uniqVertices:
                        if edge.containsPoints(testVertex):
                            if i == len(surface.vertices) - 1:
                                inext = 0
                            else:
                                inext = i + 1
                            surface.vertices.insert(inext, testVertex)
                            insertedVertex = True
                            break
                    # Break out of the loop on vertices/edges too, start again at while loop
                    if insertedVertex:
                        break
        return updZonePoly

    def isEnclosedVolume(self) -> Tuple[bool, List[Surface3dEge]]:
        """Checks if the Polyhedron is enclosed, that is all its edges are used exactly twice."""
        edgeNot2orig, _ = Polyhedron.edgesNotTwoForEnclosedVolumeTest(zonePoly=self)
        if not edgeNot2orig:
            return True, []

        print("Updating Polyhedron with collinear vertices on lines")
        updatedZonePoly = self.updateZonePolygonsForMissingColinearPoints()
        edgeNot2again, _ = Polyhedron.edgesNotTwoForEnclosedVolumeTest(updatedZonePoly)
        if not edgeNot2again:
            return True, []

        return False, edgesInBoth(edgeNot2orig, edgeNot2again)

    def calcPolyhedronVolume(self) -> float:
        """Calculates the Volume of an ENCLOSED Polyhedron."""
        volume = 0.0
        for surface in self.surfaces:
            vertices = surface.vertices
            p3FaceOrigin = vertices[1] - Vertex(0, 0, 0)
            newellAreaVector = getNewellVector(vertices)
            pyramidVolume = newellAreaVector.dot(p3FaceOrigin)
            volume += pyramidVolume
        # Our newellArea vector has twice the length
        volume /= 6.0
        return volume

    def to_os_cpp_code(self):
        """For my own convenience when writting OpenStudio tests."""
        for i, sf in enumerate(self.surfaces):
            if sf.name is not None:
                name = sf.name
            else:
                name = f"Surface {i+1}"
            if name[1] == '-':
                cleaned_name = name[2:].lower().replace('-', '')
            else:
                cleaned_name = name.lower().replace('-', '')
            s = "{"
            if sf.name is not None:
                s += "\n "
            n_vertices = len(sf.vertices)
            imax = n_vertices - 1
            for i, v in enumerate(sf.vertices):
                if i > 0:
                    s += " "
                s += f"{{{v.x:+.1f}, {v.y:+.1f}, {v.z:+.1f}}}"
                if i < imax:
                    s += ",\n"
            s += "}"

            print(f'Surface {cleaned_name}({s}, m);')
            print(f'{cleaned_name}.setName("{name}");')
            print(f'{cleaned_name}.setSpace(s);\n')

    def to_eplus_cpp_code(self):
        """For my own convenience when writting EnergyPlus tests."""
        n_surfaces = len(self.surfaces)
        print(
            f"""
            Array1D_bool enteredCeilingHeight;
            state->dataGlobal->NumOfZones = 1;
            enteredCeilingHeight.dimension(state->dataGlobal->NumOfZones, false);
            state->dataHeatBal->Zone.allocate(state->dataGlobal->NumOfZones);
            state->dataHeatBal->Zone(1).HasFloor = true;
            state->dataHeatBal->Zone(1).HTSurfaceFirst = 1;
            state->dataHeatBal->Zone(1).AllSurfaceFirst = 1;
            state->dataHeatBal->Zone(1).AllSurfaceLast = {n_surfaces};

            state->dataSurface->Surface.allocate({n_surfaces});
            """
        )
        for i, sf in enumerate(self.surfaces):
            n_vertices = len(sf.vertices)
            name = sf.name
            if 'ROOF' in name:
                tilt = 0.0
                c = 'SurfaceClass::Roof'
            elif 'FLOOR' in name:
                tilt = 180.0
                c = 'SurfaceClass::Floor'
            else:
                tilt = 90.0
                c = 'SurfaceClass::Wall'

            print(
                f'''
            state->dataSurface->Surface({i+1}).Name = "{name}";
            state->dataSurface->Surface({i+1}).Sides = {n_vertices};
            state->dataSurface->Surface({i+1}).Vertex.dimension({n_vertices});
            state->dataSurface->Surface({i+1}).Class = {c};
            state->dataSurface->Surface({i+1}).Tilt = {tilt};'''
            )

            for j, v in enumerate(sf.vertices):
                print(f"    state->dataSurface->Surface({i+1}).Vertex(1) = Vector({v.x}, {v.y}, {v.z});")

__init__(surfaces)

Constructor from a list of Surface objects.

Source code in geomeffibem/polyhedron.py 
20
21
22
23
24
25
26
27
28
def __init__(self, surfaces: List[Surface]):
    """Constructor from a list of Surface objects."""
    if not isinstance(surfaces, np.ndarray) and not isinstance(surfaces, list):
        raise ValueError("Expected a list or numpy array of Surfaces")

    for i, surface in enumerate(surfaces):
        if not isinstance(surface, Surface):
            raise ValueError(f"Element {i} is not a Surface object")
    self.surfaces = surfaces

calcPolyhedronVolume()

Calculates the Volume of an ENCLOSED Polyhedron.

Source code in geomeffibem/polyhedron.py 
119
120
121
122
123
124
125
126
127
128
129
130
def calcPolyhedronVolume(self) -> float:
    """Calculates the Volume of an ENCLOSED Polyhedron."""
    volume = 0.0
    for surface in self.surfaces:
        vertices = surface.vertices
        p3FaceOrigin = vertices[1] - Vertex(0, 0, 0)
        newellAreaVector = getNewellVector(vertices)
        pyramidVolume = newellAreaVector.dot(p3FaceOrigin)
        volume += pyramidVolume
    # Our newellArea vector has twice the length
    volume /= 6.0
    return volume

edgesNotTwoForEnclosedVolumeTest(zonePoly) staticmethod

Counts the number of times an Edge is used.

Returns the ones that isn't used twice (and the ones used twice for debugging/inspection)

Source code in geomeffibem/polyhedron.py 
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
@staticmethod
def edgesNotTwoForEnclosedVolumeTest(zonePoly: Polyhedron) -> Tuple[List[Surface3dEge], List[Surface3dEge]]:
    """Counts the number of times an Edge is used.

    Returns the ones that isn't used twice (and the ones used twice for debugging/inspection)
    """
    uniqueSurface3dEdges: List[Surface3dEge] = []

    for surface in zonePoly.surfaces:
        for edge in surface.to_Surface3dEdges():
            found = False
            for uniqEdge in uniqueSurface3dEdges:
                if uniqEdge == edge:
                    uniqEdge.allSurfaces.append(surface)
                    found = True
                    break
            if not found:
                uniqueSurface3dEdges.append(edge)

    edgesNotTwoCount = [x for x in uniqueSurface3dEdges if x.count() != 2]
    edgesTwoCount = [x for x in uniqueSurface3dEdges if x.count() == 2]
    return edgesNotTwoCount, edgesTwoCount

get_surface_by_name(name)

Locate a surface by its name.

Source code in geomeffibem/polyhedron.py 
30
31
32
33
34
def get_surface_by_name(self, name):
    """Locate a surface by its name."""
    for s in self.surfaces:
        if s.name is not None and s.name == name:
            return s

isEnclosedVolume()

Checks if the Polyhedron is enclosed, that is all its edges are used exactly twice.

Source code in geomeffibem/polyhedron.py 
105
106
107
108
109
110
111
112
113
114
115
116
117
def isEnclosedVolume(self) -> Tuple[bool, List[Surface3dEge]]:
    """Checks if the Polyhedron is enclosed, that is all its edges are used exactly twice."""
    edgeNot2orig, _ = Polyhedron.edgesNotTwoForEnclosedVolumeTest(zonePoly=self)
    if not edgeNot2orig:
        return True, []

    print("Updating Polyhedron with collinear vertices on lines")
    updatedZonePoly = self.updateZonePolygonsForMissingColinearPoints()
    edgeNot2again, _ = Polyhedron.edgesNotTwoForEnclosedVolumeTest(updatedZonePoly)
    if not edgeNot2again:
        return True, []

    return False, edgesInBoth(edgeNot2orig, edgeNot2again)

numVertices()

Counts the total number of vertices for all surfaces.

Source code in geomeffibem/polyhedron.py 
36
37
38
39
40
41
def numVertices(self):
    """Counts the total number of vertices for all surfaces."""
    count = 0
    for s in self.surfaces:
        count += len(s.vertices)
    return count

to_eplus_cpp_code()

For my own convenience when writting EnergyPlus tests.

Source code in geomeffibem/polyhedron.py 
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
def to_eplus_cpp_code(self):
    """For my own convenience when writting EnergyPlus tests."""
    n_surfaces = len(self.surfaces)
    print(
        f"""
        Array1D_bool enteredCeilingHeight;
        state->dataGlobal->NumOfZones = 1;
        enteredCeilingHeight.dimension(state->dataGlobal->NumOfZones, false);
        state->dataHeatBal->Zone.allocate(state->dataGlobal->NumOfZones);
        state->dataHeatBal->Zone(1).HasFloor = true;
        state->dataHeatBal->Zone(1).HTSurfaceFirst = 1;
        state->dataHeatBal->Zone(1).AllSurfaceFirst = 1;
        state->dataHeatBal->Zone(1).AllSurfaceLast = {n_surfaces};

        state->dataSurface->Surface.allocate({n_surfaces});
        """
    )
    for i, sf in enumerate(self.surfaces):
        n_vertices = len(sf.vertices)
        name = sf.name
        if 'ROOF' in name:
            tilt = 0.0
            c = 'SurfaceClass::Roof'
        elif 'FLOOR' in name:
            tilt = 180.0
            c = 'SurfaceClass::Floor'
        else:
            tilt = 90.0
            c = 'SurfaceClass::Wall'

        print(
            f'''
        state->dataSurface->Surface({i+1}).Name = "{name}";
        state->dataSurface->Surface({i+1}).Sides = {n_vertices};
        state->dataSurface->Surface({i+1}).Vertex.dimension({n_vertices});
        state->dataSurface->Surface({i+1}).Class = {c};
        state->dataSurface->Surface({i+1}).Tilt = {tilt};'''
        )

        for j, v in enumerate(sf.vertices):
            print(f"    state->dataSurface->Surface({i+1}).Vertex(1) = Vector({v.x}, {v.y}, {v.z});")

to_os_cpp_code()

For my own convenience when writting OpenStudio tests.

Source code in geomeffibem/polyhedron.py 
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
def to_os_cpp_code(self):
    """For my own convenience when writting OpenStudio tests."""
    for i, sf in enumerate(self.surfaces):
        if sf.name is not None:
            name = sf.name
        else:
            name = f"Surface {i+1}"
        if name[1] == '-':
            cleaned_name = name[2:].lower().replace('-', '')
        else:
            cleaned_name = name.lower().replace('-', '')
        s = "{"
        if sf.name is not None:
            s += "\n "
        n_vertices = len(sf.vertices)
        imax = n_vertices - 1
        for i, v in enumerate(sf.vertices):
            if i > 0:
                s += " "
            s += f"{{{v.x:+.1f}, {v.y:+.1f}, {v.z:+.1f}}}"
            if i < imax:
                s += ",\n"
        s += "}"

        print(f'Surface {cleaned_name}({s}, m);')
        print(f'{cleaned_name}.setName("{name}");')
        print(f'{cleaned_name}.setSpace(s);\n')

uniqueVertices()

Get a list of unique vertices (uses Vertex eq operator which has a tolerance).

Source code in geomeffibem/polyhedron.py 
43
44
45
46
47
48
49
50
51
52
53
54
55
def uniqueVertices(self) -> List[Vertex]:
    """Get a list of unique vertices (uses Vertex __eq__ operator which has a tolerance)."""
    uniqueVertices: List[Vertex] = []
    for s in self.surfaces:
        for vertex in s.vertices:
            found = False
            for unique_v in uniqueVertices:
                if unique_v == vertex:
                    found = True
                    break
            if not found:
                uniqueVertices.append(vertex)
    return uniqueVertices

updateZonePolygonsForMissingColinearPoints()

Creates a new Polyhedron with extra vertices when a point is found to be on a line segment.

Source code in geomeffibem/polyhedron.py 
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
def updateZonePolygonsForMissingColinearPoints(self) -> Polyhedron:
    """Creates a new Polyhedron with extra vertices when a point is found to be on a line segment."""
    updZonePoly = copy.deepcopy(self)

    uniqVertices = self.uniqueVertices()

    for surface in updZonePoly.surfaces:
        insertedVertex = True
        while insertedVertex:
            insertedVertex = False
            for i, edge in enumerate(surface.to_Surface3dEdges()):
                for testVertex in uniqVertices:
                    if edge.containsPoints(testVertex):
                        if i == len(surface.vertices) - 1:
                            inext = 0
                        else:
                            inext = i + 1
                        surface.vertices.insert(inext, testVertex)
                        insertedVertex = True
                        break
                # Break out of the loop on vertices/edges too, start again at while loop
                if insertedVertex:
                    break
    return updZonePoly

edgesInBoth(a, b)

Helper function.

Source code in geomeffibem/polyhedron.py 
203
204
205
206
207
208
209
210
211
def edgesInBoth(a: List[Surface3dEge], b: List[Surface3dEge]) -> List[Surface3dEge]:
    """Helper function."""
    in_both = []
    for edge_a in a:
        for edge_b in b:
            if edge_a == edge_b:
                in_both.append(edge_a)
                break
    return in_both