From e3081538f71c9eba1002bd996910e4be1f3dd498 Mon Sep 17 00:00:00 2001
From: benSepanski <ben_sepanski@alumni.baylor.edu>
Date: Wed, 24 Jun 2020 16:31:23 -0500
Subject: [PATCH] Moved useful methods from mesh_topology.py into mesh.py
---
meshmode/interop/firedrake/mesh.py | 158 +++++++++++++++++++++++++++++
1 file changed, 158 insertions(+)
create mode 100644 meshmode/interop/firedrake/mesh.py
diff --git a/meshmode/interop/firedrake/mesh.py b/meshmode/interop/firedrake/mesh.py
new file mode 100644
index 00000000..d067fb27
--- /dev/null
+++ b/meshmode/interop/firedrake/mesh.py
@@ -0,0 +1,158 @@
+__copyright__ = "Copyright (C) 2020 Benjamin Sepanski"
+
+__license__ = """
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in
+all copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+THE SOFTWARE.
+"""
+
+__doc__ = """
+.. autofunction:: get_firedrake_nodal_adjacency_group
+.. autofunction:: get_firedrake_boundary_tags
+"""
+
+from warnings import warn # noqa
+import numpy as np
+
+
+# {{{ functions to extract information from Mesh Topology
+
+def get_firedrake_nodal_adjacency_group(fdrake_mesh, cells_to_use=None):
+ """
+ Create a nodal adjacency object
+ representing the nodal adjacency of :param:`fdrake_mesh`
+
+ :param fdrake_mesh: A firedrake mesh (:class:`MeshTopology`
+ or :class:`MeshGeometry`)
+ :param cells_to_use: Either
+
+ * *None*, in which case this argument is is ignored
+ * A numpy array of firedrake cell indices, in which case
+ any cell with indices not in the array :param:`cells_to_use`
+ is ignored.
+ This induces a new order on the cells.
+ The *i*th element in the returned :class:`NodalAdjacency`
+ object corresponds to the ``cells_to_use[i]``th cell
+ in the firedrake mesh.
+
+ This feature has been used when only part of the mesh
+ needs to be converted, since firedrake has no concept
+ of a "sub-mesh".
+
+ :return: A :class:`meshmode.mesh.NodalAdjacency` instance
+ representing the nodal adjacency of :param:`fdrake_mesh`
+ """
+ mesh_topology = fdrake_mesh.topology
+ # TODO... not sure how to get around the private access
+ plex = mesh_topology._plex
+
+ # dmplex cell Start/end and vertex Start/end.
+ c_start, c_end = plex.getHeightStratum(0)
+ v_start, v_end = plex.getDepthStratum(0)
+
+ # TODO... not sure how to get around the private access
+ # This maps the dmplex index of a cell to its firedrake index
+ to_fd_id = np.vectorize(mesh_topology._cell_numbering.getOffset)(
+ np.arange(c_start, c_end, dtype=np.int32))
+
+ element_to_neighbors = {}
+ verts_checked = set() # dmplex ids of vertex checked
+
+ # If using all cells, loop over them all
+ if cells_to_use is None:
+ range_ = range(c_start, c_end)
+ # Otherwise, just the ones you're using
+ else:
+ assert isinstance(cells_to_use, np.ndarray)
+ assert np.size(cells_to_use) == np.size(np.unique(cells_to_use)), \
+ "cells_to_use must have unique values"
+ assert len(np.shape(cells_to_use)) == 1 and len(cells_to_use) > 0
+ isin = np.isin(to_fd_id, cells_to_use)
+ range_ = np.arange(c_start, c_end, 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 v_start <= vert_id < v_end 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 c_start <= other_cell_id < c_end:
+ cells.append(to_fd_id[other_cell_id - c_start])
+
+ # If only using some cells, clean out extraneous ones
+ # and relabel them to new id
+ cells = set(cells)
+ if cells_to_use is not None:
+ cells = set([cells_to_use[fd_ndx] for fd_ndx in cells
+ if fd_ndx in cells_to_use])
+
+ # mark cells as neighbors
+ for cell_one in cells:
+ element_to_neighbors.setdefault(cell_one, set())
+ element_to_neighbors[cell_one] |= cells
+
+ # Count the number of cells
+ if cells_to_use is None:
+ nelements = mesh_topology.num_cells()
+ else:
+ nelements = cells_to_use.shape[0]
+
+ # Create neighbors_starts and neighbors
+ neighbors = []
+ # FIXME : Is this the right integer type to choose?
+ neighbors_starts = np.zeros(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
+ return NodalAdjacency(neighbors_starts=neighbors_starts,
+ neighbors=neighbors)
+
+
+def get_firedrake_boundary_tags(fdrake_mesh):
+ """
+ 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`)
+
+ :param fdrake_mesh: A :mod:`firedrake` :class:`MeshTopology` or
+ :class:`MeshGeometry`
+
+ :return: A tuple of boundary tags
+ """
+ from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL
+ bdy_tags = [BTAG_ALL, BTAG_REALLY_ALL]
+
+ unique_markers = fdrake_mesh.topology.exterior_facets.unique_markers
+ if unique_markers is not None:
+ bdy_tags += list(unique_markers)
+
+ return tuple(bdy_tags)
+
+# }}}
--
GitLab