Coverage for wulffpack/core/facet.py: 98%
147 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-29 11:40 +0000
1from typing import List, Tuple, Union
2import numpy as np
3from .geometry import get_angle, get_tetrahedral_volume
6class Facet(object):
7 """
8 A `Facet` carries information about a particular facet, i.e.,
9 only one facet out of potentially several that are equivalent
10 by symmetry.
12 Parameters
13 ----------
14 normal
15 The normal to the surface in Cartesian coordintes
16 energy
17 Surface energy (in units of energy per surface area)
18 rotation
19 Symmetry operation (acting on scaled coordinates) that
20 carried an arbtriraily defined "parent facet" to this facet
21 original_grain
22 If True, this facet is part of an "original grain", i.e.,
23 in the case of decahedron/icosahedron, it is part of
24 the first created tetrahedron
25 tol
26 Numerical tolerance
27 """
29 def __init__(
30 self,
31 normal: Union[List[float], Tuple[float]],
32 energy: float,
33 symmetry: np.ndarray,
34 original_grain: bool = True,
35 tol: float = 1e-5
36 ):
37 self.normal = np.array(normal) / np.linalg.norm(normal)
38 self.original_normal = self.normal.copy()
39 self.energy = energy
40 self.symmetries = [symmetry]
41 self.vertices = []
42 self.original_grain = original_grain
43 self.original_vertices = []
44 self.tol = tol
46 def __str__(self) -> str:
47 width = 42
48 s = []
49 s += ['{s:-^{n}}'.format(s=' Facet ', n=width)]
50 flds = dict(normal='', original_normal='', energy='12.4f', original_grain='')
51 for fld, fmt in flds.items():
52 s += ['{fld:16} : {value:{fmt}}'.format(fld=fld, value=getattr(self, fld), fmt=fmt)]
53 s += [f'{"symmetries":16} :']
54 for k, v in enumerate(self.symmetries):
55 s += [f' {k:4} : {v[0]}']
56 s += [f' {"":4} {v[1]}']
57 s += [f' {"":4} {v[2]}']
58 return '\n'.join(s)
60 def _repr_html_(self) -> str:
61 s = []
62 s += ['<table border="1" class="dataframe"']
63 s += [
64 '<thead><tr><th style="text-align: left;">Property</th><th>Value</th></tr></thead>'
65 ]
66 s += ['<tbody>']
67 flds = dict(normal='', original_normal='', energy='12.4f', original_grain='')
68 for fld, fmt in flds.items():
69 s += [
70 f'<tr><td style="text-align: left">{fld:22}</td>'
71 '<td>{value:{fmt}}</td><tr>'.format(value=getattr(self, fld), fmt=fmt)
72 ]
73 s += '<tr><td style="text-align: left">symmetries</td><tr>'
74 for k, v in enumerate(self.symmetries):
75 s += [f'<tr><td>{k}</td><td>{v[0]}<br>{v[1]}<br>{v[2]}</td></tr>']
76 s += ['</tbody>']
77 s += ['</table>']
78 return ''.join(s)
80 def add_vertex(self, vertex: Union[List, tuple]):
81 """
82 Add a vertex to the list of vertices. This function checks
83 that there is not already such a vertex, to ensure that
84 no vertex is duplicated. It is also checked that the vertex
85 lies in the expected plane (as defined by the first vertex and
86 the surface normal).
88 Parameters
89 ----------
90 vertex
91 The vertex in Cartesian coordinates
92 """
93 # Need to check that the vertex is not already added
94 # (relevant if more than three planes cross at the
95 # same point, such as in a non-truncated octahedron)
96 for comp_vertex in self.vertices:
97 if np.linalg.norm(vertex - comp_vertex) < self.tol:
98 return None
99 # Check that it is in the right plane
100 if self.vertices:
101 if abs(np.dot(vertex, self.normal) - self.distance_from_origin) > self.tol:
102 raise ValueError('Vertex {} corresponding to facet with normal {} '
103 'does not lie in the same plane '
104 'as previously added vertices.'.format(vertex, self.normal))
105 self.vertices.append(vertex)
107 @property
108 def distance_from_origin(self) -> float:
109 """Returns the distance from the origin to the facet."""
110 return np.dot(self.vertices[0], self.normal)
112 def remove_redundant_vertices(self):
113 """
114 Remove any vertex which is midway between two other vertices.
115 This would not work if the polygon were concave, but that is
116 not allowed anyway.
117 """
118 to_pop = []
119 for i, vertex_i in enumerate(self.vertices):
120 for j, vertex_j in enumerate(self.vertices):
121 if j <= i:
122 continue
123 v1 = vertex_j - vertex_i
124 for k, vertex_k in enumerate(self.vertices):
125 if k <= j:
126 continue
127 v2 = vertex_k - vertex_i
128 angleish = abs(get_angle(v2, v1)) / np.pi
129 if abs(angleish - round(angleish)) > self.tol / 100:
130 continue
131 # Then we have found three vertices on a line
132 # Now we need to remove the middle one
133 v3 = vertex_k - vertex_j
134 d1 = np.linalg.norm(v1)
135 d2 = np.linalg.norm(v2)
136 d3 = np.linalg.norm(v3)
137 _, indices = zip(
138 *sorted(zip([d1, d2, d3], list(range(3)))))
139 to_pop.append([k, j, i][indices[-1]])
140 for i in reversed(sorted(set(to_pop))):
141 self.vertices.pop(i)
143 @property
144 def ordered_vertices(self) -> List[np.ndarray]:
145 """
146 Returns the vertices ordered such that they enclose a polygon.
147 The first/last vertex occurs twice.
149 It is assumed (but not checked) that the coordinates lie in
150 a plane and that they form a convex polygon.
151 """
152 # Find closest pairs by comparing angles between each triplet of
153 # vertices
154 best_pairs = []
155 for i, vertex_i in enumerate(self.vertices):
156 # Its closest neighbors are those with which the largest angle is
157 # formed
158 best_angle = 0
159 for j, vertex_j in enumerate(self.vertices):
160 if i == j:
161 continue
162 for k, vertex_k in enumerate(self.vertices):
163 if k <= j or i == k:
164 continue
165 angle = get_angle(vertex_j - vertex_i, vertex_k - vertex_i)
166 if angle > best_angle:
167 best_angle = angle
168 neighbors = (j, k)
169 for neighbor in neighbors:
170 pair = tuple(sorted([i, neighbor]))
171 if pair not in best_pairs:
172 best_pairs.append(pair)
174 # Make a new list of vertices, now in order (as implicitly
175 # defined by the pairs of nearest neighbors)
176 ordered_vertices = [self.vertices[0]]
177 current_vertex_index = 0
178 while len(ordered_vertices) < len(self.vertices):
179 for pair_index, pair in enumerate(best_pairs): 179 ↛ 178line 179 didn't jump to line 178 because the loop on line 179 didn't complete
180 if current_vertex_index in pair:
181 for new_vertex_index in pair: 181 ↛ 186line 181 didn't jump to line 186 because the loop on line 181 didn't complete
182 if new_vertex_index != current_vertex_index:
183 ordered_vertices.append(
184 self.vertices[new_vertex_index])
185 break
186 best_pairs.pop(pair_index)
187 current_vertex_index = new_vertex_index
188 break
189 ordered_vertices.append(self.vertices[0])
190 return ordered_vertices
192 @property
193 def face_as_triangles(self) -> List[List[np.ndarray]]:
194 """
195 Given N 3D coordinates, ABCD...N, split them into triangles
196 that do not overlap. Two of the vertices of each triangle
197 will be vertices of the face, the third vertex will be
198 the origin.
199 """
200 vertices = self.ordered_vertices[:-1]
201 nvertices = len(vertices)
202 center = np.mean(vertices, axis=0) # Mid point
204 # Standardize order because some applications
205 # (such as threejs) expect that
206 normal = np.cross(vertices[1] - center,
207 vertices[0] - center)
208 reverse = (np.dot(normal, self.normal) > 0)
210 # Construct all triangles
211 triangles = []
212 for i in range(0, len(vertices)):
213 if reverse:
214 triangle = [center,
215 vertices[(i + 1) % nvertices],
216 vertices[i]]
217 else:
218 triangle = [center,
219 vertices[i],
220 vertices[(i + 1) % nvertices]]
221 triangles.append(triangle)
222 assert len(triangles) == len(vertices)
223 return triangles
225 @property
226 def area(self) -> float:
227 """Returns the total area of the facet."""
228 area = 0
229 for triangle in self.face_as_triangles:
230 v1 = triangle[1] - triangle[0]
231 v2 = triangle[2] - triangle[0]
232 area += np.linalg.norm(np.cross(v1, v2))
233 return area / 2
235 @property
236 def surface_energy(self) -> float:
237 """Returns the total surface energy of the facet."""
238 return self.area * self.energy
240 def get_volume(self, origin: Union[List, tuple] = None) -> float:
241 """
242 Calculate the volume formed by the pyramid having the face as
243 its base and the specified origin as its tip.
245 Parameters
246 ----------
247 origin : list of 3 ints
248 The origin
250 Returns
251 -------
252 The total volume
253 """
254 if origin is None: 254 ↛ 257line 254 didn't jump to line 257 because the condition on line 254 was always true
255 origin = np.array([0, 0, 0])
256 else:
257 origin = np.array(origin)
259 volume = 0
260 for triangle in self.face_as_triangles:
261 volume += get_tetrahedral_volume(triangle, origin)
262 return volume
264 @property
265 def perimeter_length(self) -> float:
266 """Returns the length of the perimeter"""
267 length = 0.0
268 ordered_vertices = self.ordered_vertices
269 for i, vertex in enumerate(ordered_vertices[1:]):
270 length += np.linalg.norm(vertex - ordered_vertices[i])
271 return length