Added a mesh topology import handler

meshmode/interop/firedrake/mesh_topology.py

+from warnings import warn  # noqa
+import numpy as np
+
+from meshmode.interop import ExternalImportHandler
+
+
+# {{{ ImportHandler for firedrake's MeshTopology class
+
+class FiredrakeMeshTopologyImporter(ExternalImportHandler):
+    """
+    An Importer for :class:`firedrake.mesh.MeshTopology`.
+    Holds the topological (as opposed to geometric) information
+    about a mesh.
+    """
+
+    def __init__(self, mesh, cells_to_use=None):
+        """
+        :param mesh: An instance :mod:`firedrake` :class:`MeshTopology` or
+                   :class:`MeshGeometry`. If an instance of
+                   :class:`MeshGeometry`, uses its underlying topology.
+        :param cells_to_use: Either
+
+            * *None*, in which case this argument is ignored
+            * An array of cell ids, in which case those are the
+              only cells for which information is gathered/converted
+
+        We require that :param:`mesh` have co-dimesnion
+        of 0 or 1.
+        Moreover, if :param:`mesh` is a 2-surface embedded in 3-space,
+        we _require_ that :function:`init_cell_orientations`
+        has been called already.
+
+        :raises TypeError: If :param:`mesh` is not of :mod:`firedrake`
+                           :class:`MeshTopology` or :class:`MeshGeometry`
+        """
+        top = mesh.topological  # convert geometric to topological
+
+        # {{{ Check input types
+        from firedrake.mesh import MeshTopology
+        if not isinstance(top, MeshTopology):
+            raise TypeError(":param:`mesh` must be of type "
+                            ":class:`firedrake.mesh.MeshTopology` or "
+                            ":class:`firedrake.mesh.MeshGeometry`")
+        # }}}
+
+        super(FiredrakeMeshTopologyImporter, self).__init__(top)
+
+        # Ensure has simplex-type elements
+        if not top.ufl_cell().is_simplex():
+            raise ValueError(":param:`mesh` must have simplex type elements, "
+                             "%s is not a simplex" % (mesh.ufl_cell()))
+
+        # Ensure dimensions are in appropriate ranges
+        supported_dims = [1, 2, 3]
+        if self.cell_dimension() not in supported_dims:
+            raise ValueError("Cell dimension is %s. Cell dimension must be one of"
+                             " %s" % (self.cell_dimension(), supported_dims))
+
+        self._nodal_adjacency = None
+        self.icell_to_fd = cells_to_use  # Map cell index -> fd cell index
+        self.fd_to_icell = None          # Map fd cell index -> cell index
+
+        # perform checks on :param:`cells_to_use` if not *None*
+        if self.icell_to_fd is not None:
+            assert np.unique(self.icell_to_fd).shape == self.icell_to_fd.shape
+            self.fd_to_icell = dict(zip(self.icell_to_fd,
+                                        np.arange(self.icell_to_fd.shape[0],
+                                                  dtype=np.int32)
+                                        ))
+
+    @property
+    def topology_importer(self):
+        """
+        A reference to *self*, for compatability with mesh geometry
+        importers
+        """
+        return self
+
+    @property
+    def topological_importer(self):
+        """
+        A reference to *self*, for compatability with mesh geometry
+        importers
+        """
+        return self
+
+    def cell_dimension(self):
+        """
+        Return the dimension of the cells used by this topology
+        """
+        return self.data.cell_dimension()
+
+    def nelements(self):
+        """
+        Return the number of cells in this topology
+        """
+        if self.icell_to_fd is None:
+            num_cells = self.data.num_cells()
+        else:
+            num_cells = self.icell_to_fd.shape[0]
+
+        return num_cells
+
+    def nunit_vertices(self):
+        """
+        Return the number of unit vertices on the reference element
+        """
+        return self.data.ufl_cell().num_vertices()
+
+    def bdy_tags(self):
+        """
+        Return a tuple of bdy tags as requested in
+        the construction of a :mod:`meshmode` :class:`Mesh`
+
+        The tags used are :class:`meshmode.mesh.BTAG_ALL`,
+        :class:`meshmode.mesh.BTAG_REALLY_ALL`, and
+        any markers in the mesh topology's exterior facets
+        (see :attr:`firedrake.mesh.MeshTopology.exterior_facets.unique_markers`)
+        """
+        from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL
+        bdy_tags = [BTAG_ALL, BTAG_REALLY_ALL]
+
+        unique_markers = self.data.exterior_facets.unique_markers
+        if unique_markers is not None:
+            bdy_tags += list(unique_markers)
+
+        return tuple(bdy_tags)
+
+    def nodal_adjacency(self):
+        """
+        Returns a :class:`meshmode.mesh.NodalAdjacency` object
+        representing the nodal adjacency of this mesh
+        """
+        if self._nodal_adjacency is None:
+            # TODO... not sure how to get around the private access
+            plex = self.data._plex
+
+            # variable names follow a dmplex naming convention for
+            # cell Start/end and vertex Start/end.
+            cStart, cEnd = plex.getHeightStratum(0)
+            vStart, vEnd = plex.getDepthStratum(0)
+
+            # TODO... not sure how to get around the private access
+            to_fd_id = np.vectorize(self.data._cell_numbering.getOffset)(
+                np.arange(cStart, cEnd, dtype=np.int32))
+
+            element_to_neighbors = {}
+            verts_checked = set()  # dmplex ids of vertex checked
+
+            # If using all cells, loop over them all
+            if self.icell_to_fd is None:
+                range_ = range(cStart, cEnd)
+            # Otherwise, just the ones you're using
+            else:
+                isin = np.isin(to_fd_id, self.icell_to_fd)
+                range_ = np.arange(cStart, cEnd, dtype=np.int32)[isin]
+
+            # For each cell
+            for cell_id in range_:
+                # For each vertex touching the cell (that haven't already seen)
+                for vert_id in plex.getTransitiveClosure(cell_id)[0]:
+                    if vStart <= vert_id < vEnd and vert_id not in verts_checked:
+                        verts_checked.add(vert_id)
+                        cells = []
+                        # Record all cells touching that vertex
+                        support = plex.getTransitiveClosure(vert_id,
+                                                            useCone=False)[0]
+                        for other_cell_id in support:
+                            if cStart <= other_cell_id < cEnd:
+                                cells.append(to_fd_id[other_cell_id - cStart])
+
+                        # If only using some cells, clean out extraneous ones
+                        # and relabel them to new id
+                        cells = set(cells)
+                        if self.fd_to_icell is not None:
+                            cells = set([self.fd_to_icell[fd_ndx]
+                                         for fd_ndx in cells
+                                         if fd_ndx in self.fd_to_icell])
+
+                        # mark cells as neighbors
+                        for cell_one in cells:
+                            element_to_neighbors.setdefault(cell_one, set())
+                            element_to_neighbors[cell_one] |= cells
+
+            # Create neighbors_starts and neighbors
+            neighbors = []
+            neighbors_starts = np.zeros(self.nelements() + 1, dtype=np.int32)
+            for iel in range(len(element_to_neighbors)):
+                elt_neighbors = element_to_neighbors[iel]
+                neighbors += list(elt_neighbors)
+                neighbors_starts[iel+1] = len(neighbors)
+
+            neighbors = np.array(neighbors, dtype=np.int32)
+
+            from meshmode.mesh import NodalAdjacency
+            self._nodal_adjacency = NodalAdjacency(neighbors_starts=neighbors_starts,
+                                                   neighbors=neighbors)
+        return self._nodal_adjacency
+
+# }}}