Function build_faces

phi/geom/_mesh.py:715–831

Given a list of vertices, elements and boundary edges, computes the element connectivity matrix and corresponding edge properties. Args: vertices: `Tensor` representing list (instance) of vectors (channel) elements: Sparse matrix listing all elements (instance). Each entry

(vertices: Tensor,  # (vertices:i, vector)
                elements: Tensor,  # (elements:i, ~vertices)
                boundaries: Dict[str, Sequence],  # vertex pairs
                element_rank: int,
                periodic: Sequence[str],  # periodic dim names
                vertex_mean: Tensor,
                face_format: str)

def build_faces(vertices: Tensor,  # (vertices:i, vector)
                elements: Tensor,  # (elements:i, ~vertices)
                boundaries: Dict[str, Sequence],  # vertex pairs
                element_rank: int,
                periodic: Sequence[str],  # periodic dim names
                vertex_mean: Tensor,
                face_format: str):
    """
    Given a list of vertices, elements and boundary edges, computes the element connectivity matrix  and corresponding edge properties.

    Args:
        vertices: `Tensor` representing list (instance) of vectors (channel)
        elements: Sparse matrix listing all elements (instance). Each entry represents a vertex (dual) belonging to an element.
        boundaries: Named sequences of edges (vertex pairs).
        element_rank: Spatial rank of the elements (currently only 2 is supported)
        periodic: Which dims are periodic.
        vertex_mean: Mean vertex position for each element.
        face_format: Sparse matrix format to use for the element-element matrices.
    """
    n_v = instance(vertices).size
    n_e = instance(elements).size
    # --- Periodic: map vertices of boundary+ to the corresponding vertex in boundary- ---
    vertex_id = np.arange(instance(vertices).size)
    periodic = {dim[:-len('[::-1]')] if dim.endswith('[::-1]') else dim: dim.endswith('[::-1]') for dim in periodic}
    for dim, is_flipped in periodic.items():
        vertex_id[np.concatenate(boundaries[dim+'+'])] = vertex_id[np.concatenate(boundaries[dim+'-'])[::-1] if is_flipped else np.concatenate(boundaries[dim+'-'])]
    is_periodic = dim_mask(vertices.vector.item_names, tuple(periodic))
    # --- element-facet and facet-vertex matrix. A facet describes a single oriented face of an element, i.e. shared faces get two entries. ---
    v_count = dsum(elements).numpy()  # number of vertices per element
    v1 = stored_indices(elements).index[dual(elements).name].numpy()
    if element_rank == 2:  # edges are the lines between neighbor vertices in the vertex lists + the edge last-to-first
        v1 = vertex_id[v1]
        n_f = v1.size  # total number of facets (excluding boundaries)
        ptr = np.cumsum(v_count)
        roll = np.arange(v1.size) + 1
        roll[ptr - 1] = ptr - v_count
        v12 = np.stack([v1, v1[roll]], -1).flatten()
        f_idx = np.arange(v1.size, dtype=v1.dtype)
        f_idx2 = f_idx.repeat(2)
        f_v = coo_matrix((np.ones(n_f*2, np.int32), (f_idx2, v12)), shape=(n_f, n_v))  # facet-vertex matrix
        e_idx = np.arange(instance(elements).size).repeat(v_count)
        e_f = coo_matrix((np.ones(n_f, bool), (e_idx, f_idx)), shape=(n_e, n_f))  # element-facet matrix
        # --- Compute facet properties: center, normal, area ---
        f_v_pos = vertices[reshaped_tensor(v12, [instance('facets') + dual(pair=2)], convert=False)]  # vertex positions of every (inner) facet
        if periodic:  # map v_pos: closest to cell_center
            cell_center = vertex_mean[wrap(e_idx, 'facets:i')]
            bounds = bounding_box(vertices)
            delta = PERIODIC.shortest_distance(cell_center - bounds.lower, f_v_pos - bounds.lower, bounds.size)
            f_v_pos = where(is_periodic, cell_center + delta, f_v_pos)
        f_center = dmean(f_v_pos)
        edge_dir = f_v_pos.pair.dual[1] - f_v_pos.pair.dual[0]
        area = vec_length(edge_dir)
        normal = vec_normalize(stack([-edge_dir[1], edge_dir[0]], channel(edge_dir)))
    elif element_rank == 3:
        v3d, c3d = element_types_3d()
        n_v_per_f = [c3d[v] for v in v_count]
        n_f_per_e = [len(v) for v in n_v_per_f]
        n_fv_per_e = [sum(v) for v in n_v_per_f]