Source code for wulffpack.core.facet

from typing import List, Tuple, Union
import numpy as np
from .geometry import get_angle, get_tetrahedral_volume


[docs]class Facet(object): """ A `Facet` carries information about a particular facet, i.e., only one facet out of potentially several that are equivalent by symmetry. Parameters ---------- normal The normal to the surface in Cartesian coordintes energy Surface energy (in units of energy per surface area) rotation Symmetry operation (acting on scaled coordinates) that carried an arbtriraily defined "parent facet" to this facet original_grain If True, this facet is part of an "original grain", i.e., in the case of decahedron/icosahedron, it is part of the first created tetrahedron tol Numerical tolerance """ def __init__(self, normal: Union[List[float], Tuple[float]], energy: float, symmetry: np.ndarray, original_grain: bool=True, tol: float=1e-5): self.normal = np.array(normal) / np.linalg.norm(normal) self.original_normal = self.normal.copy() self.energy = energy self.symmetries = [symmetry] self.vertices = [] self.original_grain = original_grain self.original_vertices = [] self.tol = tol
[docs] def add_vertex(self, vertex: Union[List, tuple]): """ Add a vertex to the list of vertices. This function checks that there is not already such a vertex, to ensure that no vertex is duplicated. It is also checked that the vertex lies in the expected plane (as defined by the first vertex and the surface normal). Parameters ---------- vertex The vertex in Cartesian coordinates """ # Need to check that the vertex is not already added # (relevant if more than three planes cross at the # same point, such as in a non-truncated octahedron) for comp_vertex in self.vertices: if np.linalg.norm(vertex - comp_vertex) < self.tol: return None # Check that it is in the right plane if self.vertices: if abs(np.dot(vertex, self.normal) - self.distance_from_origin) > self.tol: raise ValueError('Vertex {} does not lie in the same plane ' 'as previously added vertices.'.format(vertex)) self.vertices.append(vertex)
@property def distance_from_origin(self) -> float: """Returns the distance from the origin to the facet.""" return np.dot(self.vertices[0], self.normal)
[docs] def remove_redundant_vertices(self): """ Remove any vertex which is midway between two other vertices. This would not work if the polygon were concave, but that is not allowed anyway. """ to_pop = [] for i, vertex_i in enumerate(self.vertices): for j, vertex_j in enumerate(self.vertices): if j <= i: continue v1 = vertex_j - vertex_i for k, vertex_k in enumerate(self.vertices): if k <= j: continue v2 = vertex_k - vertex_i angleish = abs(get_angle(v2, v1)) / np.pi if abs(angleish - round(angleish)) > self.tol / 100: continue # Then we have found three vertices on a line # Now we need to remove the middle one v3 = vertex_k - vertex_j d1 = np.linalg.norm(v1) d2 = np.linalg.norm(v2) d3 = np.linalg.norm(v3) _, indices = zip( *sorted(zip([d1, d2, d3], list(range(3))))) to_pop.append([k, j, i][indices[-1]]) for i in reversed(sorted(set(to_pop))): self.vertices.pop(i)
@property def ordered_vertices(self) -> List[np.ndarray]: """ Returns the vertices ordered such that they enclose a polygon. The first/last vertex occurs twice. It is assumed (but not checked) that the coordinates lie in a plane and that they form a convex polygon. """ # Find closest pairs by comparing angles between each triplet of # vertices best_pairs = [] for i, vertex_i in enumerate(self.vertices): # Its closest neighbors are those with which the largest angle is # formed best_angle = 0 for j, vertex_j in enumerate(self.vertices): if i == j: continue for k, vertex_k in enumerate(self.vertices): if k <= j or i == k: continue angle = get_angle(vertex_j - vertex_i, vertex_k - vertex_i) if angle > best_angle: best_angle = angle neighbors = (j, k) for neighbor in neighbors: pair = tuple(sorted([i, neighbor])) if pair not in best_pairs: best_pairs.append(pair) # Make a new list of vertices, now in order (as implicitly # defined by the pairs of nearest neighbors) ordered_vertices = [self.vertices[0]] current_vertex_index = 0 while len(ordered_vertices) < len(self.vertices): for pair_index, pair in enumerate(best_pairs): if current_vertex_index in pair: for new_vertex_index in pair: if new_vertex_index != current_vertex_index: ordered_vertices.append( self.vertices[new_vertex_index]) break best_pairs.pop(pair_index) current_vertex_index = new_vertex_index break ordered_vertices.append(self.vertices[0]) return ordered_vertices @property def face_as_triangles(self) -> List[List[np.ndarray]]: """ Given N 3D coordinates, ABCD...N, split them into triangles that do not overlap. Two of the vertices of each triangle will be vertices of the face, the third vertex will be the origin. """ vertices = self.ordered_vertices[:-1] nvertices = len(vertices) center = np.mean(vertices, axis=0) # Mid point # Standardize order because some applications # (such as threejs) expect that normal = np.cross(vertices[1] - center, vertices[0] - center) reverse = (np.dot(normal, self.normal) > 0) # Construct all triangles triangles = [] for i in range(0, len(vertices)): if reverse: triangle = [center, vertices[(i + 1) % nvertices], vertices[i]] else: triangle = [center, vertices[i], vertices[(i + 1) % nvertices]] triangles.append(triangle) assert len(triangles) == len(vertices) return triangles @property def area(self) -> float: """Returns the total area of the facet.""" area = 0 for triangle in self.face_as_triangles: v1 = triangle[1] - triangle[0] v2 = triangle[2] - triangle[0] area += np.linalg.norm(np.cross(v1, v2)) return area / 2 @property def surface_energy(self) -> float: """Returns the total surface energy of the facet.""" return self.area * self.energy
[docs] def get_volume(self, origin: Union[List, tuple]=None) -> float: """ Calculate the volume formed by the pyramid having the face as its base and the specified origin as its tip. Parameters ---------- origin : list of 3 ints The origin Returns ------- The total volume """ if origin is None: origin = np.array([0, 0, 0]) else: origin = np.array(origin) volume = 0 for triangle in self.face_as_triangles: volume += get_tetrahedral_volume(triangle, origin) return volume
@property def perimeter_length(self) -> float: """Returns the length of the perimeter""" length = 0.0 ordered_vertices = self.ordered_vertices for i, vertex in enumerate(ordered_vertices[1:]): length += np.linalg.norm(vertex - ordered_vertices[i]) return length