Coverage for wulffpack/core/facet.py: 98%
120 statements
« prev ^ index » next coverage.py v7.4.4, created at 2024-04-02 08:51 +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 add_vertex(self, vertex: Union[List, tuple]):
47 """
48 Add a vertex to the list of vertices. This function checks
49 that there is not already such a vertex, to ensure that
50 no vertex is duplicated. It is also checked that the vertex
51 lies in the expected plane (as defined by the first vertex and
52 the surface normal).
54 Parameters
55 ----------
56 vertex
57 The vertex in Cartesian coordinates
58 """
59 # Need to check that the vertex is not already added
60 # (relevant if more than three planes cross at the
61 # same point, such as in a non-truncated octahedron)
62 for comp_vertex in self.vertices:
63 if np.linalg.norm(vertex - comp_vertex) < self.tol:
64 return None
65 # Check that it is in the right plane
66 if self.vertices:
67 if abs(np.dot(vertex, self.normal) - self.distance_from_origin) > self.tol:
68 raise ValueError('Vertex {} corresponding to facet with normal {} '
69 'does not lie in the same plane '
70 'as previously added vertices.'.format(vertex, self.normal))
71 self.vertices.append(vertex)
73 @property
74 def distance_from_origin(self) -> float:
75 """Returns the distance from the origin to the facet."""
76 return np.dot(self.vertices[0], self.normal)
78 def remove_redundant_vertices(self):
79 """
80 Remove any vertex which is midway between two other vertices.
81 This would not work if the polygon were concave, but that is
82 not allowed anyway.
83 """
84 to_pop = []
85 for i, vertex_i in enumerate(self.vertices):
86 for j, vertex_j in enumerate(self.vertices):
87 if j <= i:
88 continue
89 v1 = vertex_j - vertex_i
90 for k, vertex_k in enumerate(self.vertices):
91 if k <= j:
92 continue
93 v2 = vertex_k - vertex_i
94 angleish = abs(get_angle(v2, v1)) / np.pi
95 if abs(angleish - round(angleish)) > self.tol / 100:
96 continue
97 # Then we have found three vertices on a line
98 # Now we need to remove the middle one
99 v3 = vertex_k - vertex_j
100 d1 = np.linalg.norm(v1)
101 d2 = np.linalg.norm(v2)
102 d3 = np.linalg.norm(v3)
103 _, indices = zip(
104 *sorted(zip([d1, d2, d3], list(range(3)))))
105 to_pop.append([k, j, i][indices[-1]])
106 for i in reversed(sorted(set(to_pop))):
107 self.vertices.pop(i)
109 @property
110 def ordered_vertices(self) -> List[np.ndarray]:
111 """
112 Returns the vertices ordered such that they enclose a polygon.
113 The first/last vertex occurs twice.
115 It is assumed (but not checked) that the coordinates lie in
116 a plane and that they form a convex polygon.
117 """
118 # Find closest pairs by comparing angles between each triplet of
119 # vertices
120 best_pairs = []
121 for i, vertex_i in enumerate(self.vertices):
122 # Its closest neighbors are those with which the largest angle is
123 # formed
124 best_angle = 0
125 for j, vertex_j in enumerate(self.vertices):
126 if i == j:
127 continue
128 for k, vertex_k in enumerate(self.vertices):
129 if k <= j or i == k:
130 continue
131 angle = get_angle(vertex_j - vertex_i, vertex_k - vertex_i)
132 if angle > best_angle:
133 best_angle = angle
134 neighbors = (j, k)
135 for neighbor in neighbors:
136 pair = tuple(sorted([i, neighbor]))
137 if pair not in best_pairs:
138 best_pairs.append(pair)
140 # Make a new list of vertices, now in order (as implicitly
141 # defined by the pairs of nearest neighbors)
142 ordered_vertices = [self.vertices[0]]
143 current_vertex_index = 0
144 while len(ordered_vertices) < len(self.vertices):
145 for pair_index, pair in enumerate(best_pairs): 145 ↛ 144line 145 didn't jump to line 144, because the loop on line 145 didn't complete
146 if current_vertex_index in pair:
147 for new_vertex_index in pair: 147 ↛ 152line 147 didn't jump to line 152, because the loop on line 147 didn't complete
148 if new_vertex_index != current_vertex_index:
149 ordered_vertices.append(
150 self.vertices[new_vertex_index])
151 break
152 best_pairs.pop(pair_index)
153 current_vertex_index = new_vertex_index
154 break
155 ordered_vertices.append(self.vertices[0])
156 return ordered_vertices
158 @property
159 def face_as_triangles(self) -> List[List[np.ndarray]]:
160 """
161 Given N 3D coordinates, ABCD...N, split them into triangles
162 that do not overlap. Two of the vertices of each triangle
163 will be vertices of the face, the third vertex will be
164 the origin.
165 """
166 vertices = self.ordered_vertices[:-1]
167 nvertices = len(vertices)
168 center = np.mean(vertices, axis=0) # Mid point
170 # Standardize order because some applications
171 # (such as threejs) expect that
172 normal = np.cross(vertices[1] - center,
173 vertices[0] - center)
174 reverse = (np.dot(normal, self.normal) > 0)
176 # Construct all triangles
177 triangles = []
178 for i in range(0, len(vertices)):
179 if reverse:
180 triangle = [center,
181 vertices[(i + 1) % nvertices],
182 vertices[i]]
183 else:
184 triangle = [center,
185 vertices[i],
186 vertices[(i + 1) % nvertices]]
187 triangles.append(triangle)
188 assert len(triangles) == len(vertices)
189 return triangles
191 @property
192 def area(self) -> float:
193 """Returns the total area of the facet."""
194 area = 0
195 for triangle in self.face_as_triangles:
196 v1 = triangle[1] - triangle[0]
197 v2 = triangle[2] - triangle[0]
198 area += np.linalg.norm(np.cross(v1, v2))
199 return area / 2
201 @property
202 def surface_energy(self) -> float:
203 """Returns the total surface energy of the facet."""
204 return self.area * self.energy
206 def get_volume(self, origin: Union[List, tuple] = None) -> float:
207 """
208 Calculate the volume formed by the pyramid having the face as
209 its base and the specified origin as its tip.
211 Parameters
212 ----------
213 origin : list of 3 ints
214 The origin
216 Returns
217 -------
218 The total volume
219 """
220 if origin is None: 220 ↛ 223line 220 didn't jump to line 223, because the condition on line 220 was never false
221 origin = np.array([0, 0, 0])
222 else:
223 origin = np.array(origin)
225 volume = 0
226 for triangle in self.face_as_triangles:
227 volume += get_tetrahedral_volume(triangle, origin)
228 return volume
230 @property
231 def perimeter_length(self) -> float:
232 """Returns the length of the perimeter"""
233 length = 0.0
234 ordered_vertices = self.ordered_vertices
235 for i, vertex in enumerate(ordered_vertices[1:]):
236 length += np.linalg.norm(vertex - ordered_vertices[i])
237 return length