Merge pull request #196 from tataratat/ft-expand

Ft expand

gustaf/_version.py

@@ -3,4 +3,4 @@
 Current version.
 """
 
-version = "0.0.23"
+version = "0.0.24"

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,149 @@ def to_quad(tri):
     faces[:, 3] = edge_mid_column.reshape(-1, 3)[:, [2, 0, 1]].ravel()
 
     return Faces(vertices=vertices, faces=faces)
+
+
+def edges_as_quad(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 vertex_normals(
+    faces,
+    area_weighting=False,
+    angle_weighting=False,
+    return_original_ids=False,
+):
+    """
+    Computes vertex normals and saves it in vertex_data.
+    This calls inplace remove_unreferenced_vertices, but original IDs can be
+    retrieved using the flag `return_original_ids`.
+
+    The normals are computed on the face-centers and their contributions are
+    weighted and added to the vertex normals. Per default, all element faces
+    that are adjacent to a vertex are added with equal contributions, but it is
+    also possible to use weightings by area of the adjacent element
+    (`area_weighting`) or by the angle between edges at the corner vertex.
+
+    Parameters
+    ----------
+    faces: Faces
+    area_weighting : bool (false)
+      Use the element area as a weighting to its respective normal contribution
+    angle_weighting : bool (false)
+      Use the angle of between element edges as a weighting to its respective
+      normal contribution
+    return_original_ids : bool (false)
+      return the original ids in the global mesh
+
+    Returns
+    -------
+    faces: Faces
+      faces with vertex_data["normals"] computed.
+    """
+    if not faces.whatami.startswith("tri"):
+        raise ValueError("Vertex normals only supports triangle faces")
+
+    if faces.vertices.shape[1] != 3:
+        raise ValueError("Vertex normals only support 3d triangles")
+
+    if return_original_ids:
+        original_ids = np.where(faces.referenced_vertices())[0]
+
+    faces.remove_unreferenced_vertices()
+
+    triangles = faces.vertices[faces.faces]
+
+    # compute (1 - 0) and (2 - 0)
+    edge_ab = triangles[:, 1] - triangles[:, 0]
+    edge_bc = triangles[:, 2] - triangles[:, 1]
+
+    # get normal of each faces and normalize
+    crossed = utils.arr.cross3d(edge_ab, edge_bc)
+    crossed_length = np.linalg.norm(crossed, axis=1).reshape(-1, 1)
+
+    weights = np.ones_like(crossed, dtype=settings.FLOAT_DTYPE)
+
+    # get area based weights
+    if area_weighting:
+        weights *= crossed_length
+    #     crossed /= crossed_length
+    # else:
+    crossed /= crossed_length
+
+    # get triangle corner angles (same as faces.faces)
+    if angle_weighting:
+        angles = np.empty_like(crossed, dtype=settings.FLOAT_DTYPE)
+        norm_ab = np.linalg.norm(edge_ab, axis=1)
+        norm_bc = np.linalg.norm(edge_bc, axis=1)
+        edge_ca = edge_ab + edge_bc
+        norm_ca = np.linalg.norm(edge_ca, axis=1)
+        np.arcsin(
+            crossed_length.ravel() / (norm_ab * norm_bc).ravel(),
+            out=angles[:, 0],
+        )
+        np.arccos(
+            np.einsum("ij,ij->i", edge_bc, edge_ca) / (norm_bc * norm_ca),
+            out=angles[:, 1],
+        )
+        angles[:, 2] = np.pi - angles[:, 0] - angles[:, 1]
+
+        weights *= angles
+
+    # initialize
+    normals = np.zeros_like(faces.vertices, dtype=settings.FLOAT_DTYPE)
+
+    # sum normals and weights
+    np.add.at(
+        normals, faces.faces[:, 0], crossed * weights[:, 0].reshape(-1, 1)
+    )
+    np.add.at(
+        normals, faces.faces[:, 1], crossed * weights[:, 1].reshape(-1, 1)
+    )
+    np.add.at(
+        normals, faces.faces[:, 2], crossed * weights[:, 2].reshape(-1, 1)
+    )
+
+    # normalize
+    normals /= np.linalg.norm(normals, axis=1).reshape(-1, 1)
+
+    faces.vertex_data["normals"] = normals
+    if return_original_ids:
+        return (faces, original_ids)
+    else:
+        return faces

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 = cross3d(der[::2], der[1::2])
+
+    if normalize:
+        normals /= 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