add expand

gustaf/create/edges.py

@@ -8,6 +8,33 @@
 from gustaf.edges import Edges
 
 
+def from_vertices(vertices, closed=False, continuous=True):
+    """
+    Creates Edges with given vertices. If close==True,
+    last vertices will be connected to the first one.
+    If continuous==False, it assumes that every two vertices form
+    an independent line.
+
+    Parameters
+    ----------
+    vertices: (n, d) np.ndarray or Vertices
+    close: bool
+    continuous: bool
+
+    Returns
+    -------
+    edges: Edges
+    """
+    if hasattr(vertices, "vertices"):  # noqa SIM108
+        v = vertices.vertices
+    else:
+        v = vertices
+
+    edges = Edges(v, utils.connec.range_to_edges(len(v), closed, continuous))
+
+    return edges
+
+
 def from_data(gus_obj, data, scale=None, data_norm=None):
     """
     Creates edges from gustaf object with vertices.

gustaf/create/faces.py

@@ -203,3 +203,104 @@ def to_quad(tri):
     faces[:, 3] = edge_mid_column.reshape(-1, 3)[:, [2, 0, 1]].ravel()
 
     return Faces(vertices=vertices, faces=faces)
+
+
+def expand(edges, scaled_normals):
+    """
+    Expands edges to hexa with given scaled_normals.
+
+    Parameters
+    ----------
+    edges: Edges
+      (n, d) vertices, (m, 2) edges
+    scaled_normals: (n, d) np.ndarray
+      Values will be used to subtract/add to existing vertices.
+
+    Returns
+    -------
+    expanded: Faces
+    """
+    if edges.whatami != "edges":
+        raise TypeError("expand currently only supports Edges.")
+
+    size = scaled_normals.shape[0]
+    dim = scaled_normals.shape[1]
+    if len(edges.vertices) != size:
+        raise ValueError("Edges.vertices and scaled_normals size mismatch.")
+
+    if dim != edges.vertices.shape[1]:
+        raise ValueError(
+            "Dimension mismatch between Edges.Vertices and scaled_normals"
+        )
+
+    vertices = np.empty((size * 2, dim), dtype=settings.FLOAT_DTYPE)
+    vertices[:size] = edges.vertices - scaled_normals
+    vertices[size:] = edges.vertices + scaled_normals
+
+    quad = np.empty((len(edges.edges), 4), dtype=settings.INT_DTYPE)
+    quad[:, :2] = edges.edges
+    quad[:, 2] = edges.edges[:, 1] + size
+    quad[:, 3] = edges.edges[:, 0] + size
+
+    return Faces(vertices, quad)
+
+
+def _two_ortho(trajectory):
+    aux = np.empty_like(trajectory)
+    aux[:, 0] = trajectory[:, -2]
+    aux[:, [1, 2]] = trajectory[:, [0, 1]]
+
+    # cross and normalize
+    o1 = utils.arr.cross3d(trajectory, aux)
+    norm_inv = np.linalg.norm(o1, axis=1)
+    np.reciprocal(norm_inv, out=norm_inv)
+    np.multiply(o1, norm_inv.reshape(-1, 1), out=o1)
+
+    # cross and normalize
+    o2 = np.cross(trajectory, o1)
+    norm_inv = np.linalg.norm(o2, axis=1)
+    np.reciprocal(norm_inv, out=norm_inv)
+    np.multiply(o2, norm_inv.reshape(-1, 1), out=o2)
+
+    return o1, o2
+
+
+def _circle_sample(centers, r, o1, o2, resolution):
+    """given two orthogonal vectors, this samples a circle on that plane"""
+    o_shape = o1.shape
+
+    cir = np.empty((resolution * o_shape[0], o_shape[1]))
+
+    q = np.linspace(0, 2 * np.pi, resolution + 1)[:-1]
+
+    cos = np.cos(q)
+    sin = np.sin(q)
+
+    for i, (c, s) in enumerate(zip(cos, sin)):
+        cir[i::resolution] = centers + (r * c * o1) + (r * s * o2)
+
+    return cir
+
+
+def cylinder_expand(edges, r, trajectory, resolution=10):
+    vertices = _circle_sample(
+        edges.vertices, r, *_two_ortho(trajectory), resolution
+    )
+
+    e = edges.edges * resolution
+    flip_e = e[:, [1, 0]]
+
+    quads = np.empty(
+        (len(edges.edges) * resolution, 4), dtype=settings.INT_DTYPE
+    )
+    quads[::resolution, [0, 1]] = e
+    quads[::resolution, [2, 3]] = flip_e + 1
+    for i in range(1, resolution - 1):
+        # offset_i = offset * i
+        quads[i::resolution, [0, 1]] = e + i
+        quads[i::resolution, [2, 3]] = flip_e + (i + 1)
+
+    quads[(resolution - 1) :: resolution, [0, 1]] = e + (resolution - 1)
+    quads[(resolution - 1) :: resolution, [2, 3]] = flip_e
+
+    return Faces(vertices, quads)

gustaf/utils/arr.py

@@ -525,3 +525,83 @@ def is_one_of_shapes(arr, shapes, strict=False):
         return False
 
     return True
+
+
+def derivatives_to_normals(derivatives, normalize=True):
+    """
+    Parameters
+    ----------
+    derivatives: (n, (d - 1), d) np.ndarray
+      Surface jacobian transposed.
+    normalize: bool
+
+    Returns
+    -------
+    normals: (n, d) np.ndarray
+    """
+    if derivatives.ndim != 3:
+        raise ValueError("derivatives for normals expect 3D arrays")
+
+    shape = derivatives.shape
+    if shape[0] != shape[1] - 1:
+        raise ValueError("derivatives are expected to have (d-1, d) shape")
+
+    # 2D is simple index flip
+    if shape[2] == 2:
+        der = derivatives.reshape(-1, shape[2])
+        normals = np.empty_like(der)
+        normals[:, 0] = der[:, 1]
+        normals[:, 1] = -der[:, 0]
+
+    elif shape[2] == 3:
+        der = derivatives.reshape(shape[0] * shape[1], shape[2])
+        normals = np.cross(der[::2], der[1::2])
+
+    if normalize:
+        normals *= 1.0 / np.linalg.norm(normals, axis=1).reshape(-1, 1)
+
+    return normals
+
+
+def cross3d(a, b):
+    """
+    Cross product for two 3D arrays. Usually faster than np.cross
+    as it just targets 3d.
+
+    Parameters
+    ----------
+    a: (n, 3) np.ndarray
+    b: (n, 3) np.ndarray
+
+    Returns
+    -------
+    crossed: (n, 3) np.ndarray
+    """
+    # (1 5 - 2 4, 2 3 - 0 5, 0 4 - 1 3).
+    # or from two arrays
+    # (1 2 - 2 1, 2 0 - 0 2, 0 1 - 1 0).
+    o = np.empty_like(a)
+
+    # temporary aux arrays
+    size = len(a)
+    t0 = np.empty(size)
+    t1 = np.empty(size)
+
+    # short cuts
+    a0, a1, a2 = a[..., 0], a[..., 1], a[..., 2]
+    b0, b1, b2 = b[..., 0], b[..., 1], b[..., 2]
+    o0, o1, o2 = o[..., 0], o[..., 1], o[..., 2]
+
+    np.multiply(a1, b2, out=t0)
+    np.multiply(a2, b1, out=t1)
+    np.subtract(t0, t1, out=o0)
+
+    np.multiply(a2, b0, out=t0)
+    np.multiply(a0, b2, out=t1)
+    np.subtract(t0, t1, out=o1)
+
+    np.multiply(a0, b1, out=t0)
+    np.multiply(a1, b0, out=t1)
+    np.subtract(t0, t1, out=o2)
+
+    return o