diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 4bf68cce18cf9ae4d3e72d5a8b4420b5085744aa..96a0e587fb572f5c034bc0ff38bba6650482f291 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -24,6 +24,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. """ +import six import numpy as np import modepy as mp import numpy.linalg as la @@ -40,12 +41,56 @@ __doc__ = """ :undoc-members: .. autoclass:: NodalAdjacency +.. autoclass:: FacialAdjacencyGroup .. autofunction:: as_python +Predefined Boundary tags +------------------------ + +.. autoclass:: BTAG_NONE +.. autoclass:: BTAG_ALL +.. autoclass:: BTAG_REALLY_ALL +.. autoclass:: BTAG_NO_BOUNDARY """ +# {{{ element tags + +class BTAG_NONE(object): + """A boundary tag representing an empty boundary or volume.""" + pass + + +class BTAG_ALL(object): + """A boundary tag representing the entire boundary or volume. + + In the case of the boundary, TAG_ALL does not include rank boundaries, + or, more generally, anything tagged with TAG_NO_BOUNDARY.""" + pass + + +class BTAG_REALLY_ALL(object): + """A boundary tag representing the entire boundary. + + Unlike :class:`TAG_ALL`, this includes rank boundaries, + or, more generally, everything tagged with :class:`TAG_NO_BOUNDARY`.""" + pass + + +class BTAG_NO_BOUNDARY(object): + """A boundary tag indicating that this edge should not fall under + :class:`TAG_ALL`. Among other things, this is used to keep rank boundaries + out of :class:`BTAG_ALL`. + """ + pass + + +SYSTEM_TAGS = set([BTAG_NONE, BTAG_ALL, BTAG_REALLY_ALL, BTAG_NO_BOUNDARY]) + +# }}} + + # {{{ element group class MeshElementGroup(Record): @@ -237,15 +282,15 @@ class SimplexElementGroup(MeshElementGroup): # }}} -# {{{ mesh +# {{{ nodal adjacency class NodalAdjacency(Record): - """Describes element adjacency information, i.e. information about + """Describes nodal element adjacency information, i.e. information about elements that touch in at least one point. .. attribute:: neighbors_starts - ``element_id_t [nelments+1]`` + ``element_id_t [nelements+1]`` Use together with :attr:`neighbors`. ``neighbors_starts[iel]`` and ``neighbors_starts[iel+1]`` together indicate a ranges of element indices @@ -271,6 +316,76 @@ class NodalAdjacency(Record): def __ne__(self, other): return not self.__eq__(other) +# }}} + + +# {{{ facial adjacency + +class FacialAdjacencyGroup(Record): + """Describes facial element adjacency information for one + :class:`MeshElementGroup`, i.e. information about elements that share (part + of) a face. + + .. attribute:: igroup + + .. attribute:: ineighbor_group + + ID of neighboring group, or *None* for boundary faces. If identical + to :attr:`igroup`, then this contains the self-connectivity in this + group. + + .. attribute:: elements + + ``element_id_t [nfagrp_elements]``. ``elements[i]`` + + .. attribute:: element_faces + + ``face_id_t [nfagrp_elements]``. ``element_faces[i]`` + indicate what face index of the opposite element indicated in + ``neighbors[iel_grp][iface]`` touches face number *iface* of element + number *iel_grp* in this element group. + + .. attribute:: neighbors + + ``element_id_t [nfagrp_elements]``. ``neighbors[i]`` + gives the element number within :attr:`ineighbor_group` of the element + opposite ``elements[i]``. + + If this number is negative, then this indicates that this is a + boundary face, and the bits set in ``-neighbors[i]`` + should be interpreted according to :attr:`Mesh.boundary_tags`. + + .. attribute:: neighbor_faces + + ``face_id_t [nfagrp_elements]``. ``neighbor_faces[i]`` indicate what + face index of the opposite element indicated in ``neighbors[i]`` has + facial contact with face number ``element_faces[i]`` of element number + ``elements[i]`` in element group *igroup*. + + Zero if ``neighbors[i]`` is negative. + + .. automethod:: __eq__ + .. automethod:: __ne__ + """ + + def __eq__(self, other): + return ( + type(self) == type(other) + and self.igroup == other.igroup + and self.ineighbor_group == other.ineighbor_group + and np.array_equal(self.elements, other.elements) + and np.array_equal(self.element_faces, other.element_faces) + and np.array_equal(self.neighbors, other.neighbors) + and np.array_equal(self.neighbor_faces, other.neighbor_faces) + ) + + def __ne__(self, other): + return not self.__eq__(other) + +# }}} + + +# {{{ mesh class Mesh(Record): """ @@ -290,6 +405,30 @@ class Mesh(Record): Referencing this attribute may raise :exc:`meshmode.DataUnavailable`. + .. attribute:: facial_adjacency_groups + + An list of mappings from neighbor group IDs to instances of + :class:`FacialAdjacencyGroup`. + + ``facial_adjacency_groups[igrp][ineighbor_group]`` gives + the set of facial adjacency relations between group *igrp* + and *ineighbor_group*. *ineighbor_group* and *igrp* may be + identical, or *ineighbor_group* may be *None*, in which case + a group containing boundary faces is returned. + + Referencing this attribute may raise + :exc:`meshmode.DataUnavailable`. + + .. attribute:: boundary_tags + + A tuple of boundary tag identifiers. :class:`BTAG_ALL` and + :class:`BTAG_REALLY_ALL` are guranateed to exist. + + .. attribute:: btag_to_index + + A mapping that maps boundary tag identifiers to their + corresponding index. + .. attribute:: vertex_id_dtype .. attribute:: element_id_dtype @@ -298,8 +437,12 @@ class Mesh(Record): .. automethod:: __ne__ """ + face_id_dtype = np.int8 + def __init__(self, vertices, groups, skip_tests=False, nodal_adjacency=False, + facial_adjacency_groups=False, + boundary_tags=None, vertex_id_dtype=np.int32, element_id_dtype=np.int32): """ @@ -315,8 +458,16 @@ class Mesh(Record): the full picture), and references to :attr:`element_neighbors_starts` and :attr:`element_neighbors` will result in exceptions. Lastly, a tuple - *(element_neighbors_starts, element_neighbors)*, representing the - correspondingly-named attributes. + :class:`NodalAdjacency` object. + :arg facial_adjacency_groups: One of three options: + *None*, in which case this information + will be deduced from vertex adjacency. *False*, in which case + this information will be marked unavailable (such as if there are + hanging nodes in the geometry, so that vertex adjacency does not convey + the full picture), and references to + :attr:`element_neighbors_starts` and :attr:`element_neighbors` + will result in exceptions. Lastly, a data structure as described in + :attr:`facial_adjacency_groups` may be passed. """ el_nr = 0 node_nr = 0 @@ -328,18 +479,45 @@ class Mesh(Record): el_nr += ng.nelements node_nr += ng.nnodes + # {{{ boundary tags + + if boundary_tags is None: + boundary_tags = [] + else: + boundary_tags = boundary_tags[:] + + if BTAG_ALL not in boundary_tags: + boundary_tags.append(BTAG_ALL) + if BTAG_REALLY_ALL not in boundary_tags: + boundary_tags.append(BTAG_REALLY_ALL) + + max_boundary_tag_count = int( + np.log(np.iinfo(element_id_dtype).max)/np.log(2)) + if len(boundary_tags) > max_boundary_tag_count: + raise ValueError("too few bits in element_id_dtype to represent all " + "boundary tags") + + btag_to_index = dict( + (btag, i) for i, btag in enumerate(boundary_tags)) + + # }}} + if nodal_adjacency is not False and nodal_adjacency is not None: - nb_starts, nbs = nodal_adjacency - nodal_adjacency = NodalAdjacency( - neighbors_starts=nb_starts, - neighbors=nbs) + if not isinstance(nodal_adjacency, NodalAdjacency): + nb_starts, nbs = nodal_adjacency + nodal_adjacency = NodalAdjacency( + neighbors_starts=nb_starts, + neighbors=nbs) - del nb_starts - del nbs + del nb_starts + del nbs Record.__init__( self, vertices=vertices, groups=new_groups, _nodal_adjacency=nodal_adjacency, + _facial_adjacency_groups=facial_adjacency_groups, + boundary_tags=boundary_tags, + btag_to_index=btag_to_index, vertex_id_dtype=np.dtype(vertex_id_dtype), element_id_dtype=np.dtype(element_id_dtype), ) @@ -349,6 +527,32 @@ class Mesh(Record): for g in self.groups: assert g.vertex_indices.dtype == self.vertex_id_dtype + if nodal_adjacency: + assert nodal_adjacency.neighbor_starts.shape == (self.nelements+1,) + assert len(nodal_adjacency.neighbors.shape) == 1 + + assert nodal_adjacency.neighbor_starts.dtype == self.element_id_dtype + assert nodal_adjacency.neighbors.dtype == self.element_id_dtype + + if facial_adjacency_groups: + assert len(facial_adjacency_groups) == len(self.groups) + for fagrp_map in facial_adjacency_groups: + for fagrp in six.itervalues(fagrp_map): + grp = self.groups[fagrp.igroup] + + fvi = grp.face_vertex_indices() + assert fagrp.neighbors.dtype == self.element_id_dtype + assert fagrp.neighbors.shape == ( + grp.nelements, len(fvi)) + assert fagrp.neighbors.dtype == self.face_id_dtype + assert fagrp.neighbor_faces.shape == ( + grp.nelements, len(fvi)) + + is_bdry = fagrp.neighbors < 0 + assert ((1 << btag_to_index[BTAG_REALLY_ALL]) + & fagrp.neighbors[is_bdry]).all(), \ + "boundary faces without BTAG_REALLY_ALL found" + from meshmode.mesh.processing import \ test_volume_mesh_element_orientations @@ -385,11 +589,18 @@ class Mesh(Record): passed to a Mesh constructor. """ - if isinstance(self._nodal_adjacency, NodalAdjacency): - return (self._nodal_adjacency.neighbors_starts, - self._nodal_adjacency.neighbors) - else: - return self._nodal_adjacency + return self._nodal_adjacency + + @property + def facial_adjacency_groups(self): + if self._facial_adjacency_groups is False: + from meshmode import DataUnavailable + raise DataUnavailable("facial_adjacency_groups") + elif self._facial_adjacency_groups is None: + self._facial_adjacency_groups = \ + _compute_facial_adjacency_from_vertices(self) + + return self._facial_adjacency_groups def __eq__(self, other): return ( @@ -399,7 +610,10 @@ class Mesh(Record): and self.vertex_id_dtype == other.vertex_id_dtype and self.element_id_dtype == other.element_id_dtype and (self._nodal_adjacency - == other._nodal_adjacency)) + == other._nodal_adjacency) + and (self._facial_adjacency_groups + == other._facial_adjacency_groups) + and self.boundary_tags == other.boundary_tags) def __ne__(self, other): return not self.__eq__(other) @@ -461,7 +675,7 @@ def _test_node_vertex_consistency(mesh): # }}} -# {{{ vertex-based adjacency +# {{{ vertex-based nodal adjacency def _compute_nodal_adjacency_from_vertices(mesh): # FIXME Native code would make this faster @@ -503,6 +717,115 @@ def _compute_nodal_adjacency_from_vertices(mesh): # }}} +# {{{ vertex-based facial adjacency + +def _compute_facial_adjacency_from_vertices(mesh): + # FIXME Native code would make this faster + + # create face_map, which is a mapping of + # (vertices on a face) -> + # [(igrp, iel_grp, face_idx) for elements bordering that face] + face_map = {} + for igrp, grp in enumerate(mesh.groups): + for fid, face_vertex_indices in enumerate(grp.face_vertex_indices()): + all_fvi = grp.vertex_indices[:, face_vertex_indices] + + for iel_grp, fvi in enumerate(all_fvi): + face_map.setdefault( + frozenset(fvi), []).append((igrp, iel_grp, fid)) + + # maps tuples (igrp, ineighbor_group) to number of elements + group_count = {} + for face_tuples in six.itervalues(face_map): + if len(face_tuples) == 2: + (igrp, _, _), (inb_grp, _, _) = face_tuples + group_count[igrp, inb_grp] = group_count.get((igrp, inb_grp), 0) + 1 + group_count[inb_grp, igrp] = group_count.get((inb_grp, igrp), 0) + 1 + elif len(face_tuples) == 1: + (igrp, _, _), = face_tuples + group_count[igrp, None] = group_count.get((igrp, None), 0) + 1 + else: + raise RuntimeError("unexpected number of adjacent faces") + + # {{{ build facial_adjacency_groups data structure, still empty + + facial_adjacency_groups = [] + for igroup in range(len(mesh.groups)): + grp_map = {} + facial_adjacency_groups.append(grp_map) + + bdry_count = group_count.get((igroup, None)) + if bdry_count is not None: + elements = np.empty(bdry_count, dtype=mesh.element_id_dtype) + element_faces = np.empty(bdry_count, dtype=mesh.face_id_dtype) + neighbors = np.empty(bdry_count, dtype=mesh.element_id_dtype) + neighbor_faces = np.zeros(bdry_count, dtype=mesh.face_id_dtype) + + neighbors.fill(-( + 1 << mesh.btag_to_index[BTAG_ALL] + | 1 << mesh.btag_to_index[BTAG_REALLY_ALL])) + + grp_map[None] = FacialAdjacencyGroup( + igroup=igroup, + ineighbor_group=None, + elements=elements, + element_faces=element_faces, + neighbors=neighbors, + neighbor_faces=neighbor_faces) + + for ineighbor_group in range(len(mesh.groups)): + nb_count = group_count.get((igroup, ineighbor_group)) + if nb_count is not None: + elements = np.empty(nb_count, dtype=mesh.element_id_dtype) + element_faces = np.empty(nb_count, dtype=mesh.face_id_dtype) + neighbors = np.empty(nb_count, dtype=mesh.element_id_dtype) + neighbor_faces = np.empty(nb_count, dtype=mesh.face_id_dtype) + + grp_map[ineighbor_group] = FacialAdjacencyGroup( + igroup=igroup, + ineighbor_group=ineighbor_group, + elements=elements, + element_faces=element_faces, + neighbors=neighbors, + neighbor_faces=neighbor_faces) + + # }}} + + # maps tuples (igrp, ineighbor_group) to number of elements filled in group + fill_count = {} + for face_tuples in six.itervalues(face_map): + if len(face_tuples) == 2: + for (igroup, iel, iface), (inb_group, inb_el, inb_face) in [ + (face_tuples[0], face_tuples[1]), + (face_tuples[1], face_tuples[0]), + ]: + idx = fill_count.get((igrp, inb_grp), 0) + fill_count[igrp, inb_grp] = idx + 1 + + fagrp = facial_adjacency_groups[igroup][inb_grp] + fagrp.elements[idx] = iel + fagrp.element_faces[idx] = iface + fagrp.neighbors[idx] = inb_el + fagrp.neighbor_faces[idx] = inb_face + + elif len(face_tuples) == 1: + (igroup, iel, iface), = face_tuples + + idx = fill_count.get((igrp, None), 0) + fill_count[igrp, None] = idx + 1 + + fagrp = facial_adjacency_groups[igroup][inb_grp] + fagrp.elements[idx] = iel + fagrp.element_faces[idx] = iface + + else: + raise RuntimeError("unexpected number of adjacent faces") + + return facial_adjacency_groups + +# }}} + + # {{{ as_python def _numpy_array_as_python(array): @@ -518,11 +841,15 @@ def as_python(mesh, function_name="make_mesh"): from pytools.py_codegen import PythonCodeGenerator, Indentation cg = PythonCodeGenerator() - cg("# generated by meshmode.mesh.as_python") - cg("") - cg("import numpy as np") - cg("from meshmode.mesh import Mesh, MeshElementGroup") - cg("") + cg(""" + # generated by meshmode.mesh.as_python + + import numpy as np + from meshmode.mesh import ( + Mesh, MeshElementGroup, FacialAdjacencyGroup, + BTAG_ALL, BTAG_REALLY_ALL) + + """) cg("def %s():" % function_name) with Indentation(cg): @@ -543,6 +870,46 @@ def as_python(mesh, function_name="make_mesh"): cg(" unit_nodes=%s))" % _numpy_array_as_python(group.unit_nodes)) + # {{{ facial adjacency groups + + def fagrp_params_str(fagrp): + params = { + "igroup": fagrp.igroup, + "ineighbor_group": repr(fagrp.ineighbor_group), + "elements": _numpy_array_as_python(fagrp.elements), + "element_faces": _numpy_array_as_python(fagrp.element_faces), + "neighbors": _numpy_array_as_python(fagrp.neighbors), + "neighbor_faces": _numpy_array_as_python(fagrp.neighbor_faces), + } + return ",\n ".join("%s=%s" % (k, v) for k, v in six.iteritems(params)) + + if mesh._facial_adjacency_groups: + cg("facial_adjacency_groups = []") + + for igrp, fagrp_map in enumerate(mesh.facial_adjacency_groups): + cg("facial_adjacency_groups.append({%s})" % ",\n ".join( + "%r: FacialAdjacencyGroup(%s)" % ( + inb_grp, fagrp_params_str(fagrp)) + for inb_grp, fagrp in six.iteritems(fagrp_map))) + + else: + cg("facial_adjacency_groups = %r" % mesh._facial_adjacency_groups) + + # }}} + + # {{{ boundary tags + + def strify_boundary_tag(btag): + if isinstance(btag, type): + return btag.__name__ + else: + return repr(btag) + + btags_str = ", ".join( + strify_boundary_tag(btag) for btag in mesh.boundary_tags) + + # }}} + cg("return Mesh(vertices, groups, skip_tests=True,") cg(" vertex_id_dtype=np.%s," % mesh.vertex_id_dtype.name) cg(" element_id_dtype=np.%s," % mesh.element_id_dtype.name) @@ -557,7 +924,11 @@ def as_python(mesh, function_name="make_mesh"): else: el_con_str = repr(mesh._nodal_adjacency) - cg(" nodal_adjacency=%s)" % el_con_str) + cg(" nodal_adjacency=%s," % el_con_str) + cg(" facial_adjacency_groups=facial_adjacency_groups,") + cg(" boundary_tags=[%s])" % btags_str) + + # FIXME: Handle facial adjacency, boundary tags return cg.get() diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 0f130a5d498bcdec06bb75ef0011ce9828bac6b1..61db3772bb82c41a91b1ef47221b082ba5602748 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -212,8 +212,10 @@ def make_curve_mesh(curve_f, element_boundaries, order): nodes=nodes, unit_nodes=unodes) - return Mesh(vertices=vertices, groups=[egroup], - nodal_adjacency=None) + return Mesh( + vertices=vertices, groups=[egroup], + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} @@ -285,8 +287,10 @@ def generate_icosahedron(r, order): grp = make_group_from_vertices(vertices, vertex_indices, order) from meshmode.mesh import Mesh - return Mesh(vertices, [grp], - nodal_adjacency=None) + return Mesh( + vertices, [grp], + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} @@ -302,8 +306,10 @@ def generate_icosphere(r, order): nodes=grp.nodes * r / np.sqrt(np.sum(grp.nodes**2, axis=0))) from meshmode.mesh import Mesh - return Mesh(mesh.vertices, [grp], - nodal_adjacency=None) + return Mesh( + mesh.vertices, [grp], + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} @@ -365,7 +371,11 @@ def generate_torus_and_cycle_vertices(r_outer, r_inner, nodes[2] = b*np.sin(minor_theta) from meshmode.mesh import Mesh - return (Mesh(vertices, [grp.copy(nodes=nodes)], nodal_adjacency=None), + return ( + Mesh( + vertices, [grp.copy(nodes=nodes)], + nodal_adjacency=None, + facial_adjacency_groups=None), [idx(i, 0) for i in range(n_outer)], [idx(0, j) for j in range(n_inner)]) @@ -470,7 +480,8 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64): from meshmode.mesh import Mesh return Mesh(vertices, [grp], - nodal_adjacency=None) + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 24bcae5855d29e63dc34e4b909bc660dcf29c2b5..38617cca6a6e8d6f3a215bb38c6786630783b870 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -171,7 +171,10 @@ class GmshMeshReceiver(GmshMeshReceiverBase): groups.append(group) - return Mesh(vertices, groups, nodal_adjacency=None) + return Mesh( + vertices, groups, + nodal_adjacency=None, + facial_adjacency_groups=None) # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index cfed94b0e2ef31af23977414c947cd7cf3f455fa..2bd67300e8add4be8adc1b04ff6526d3239a663b 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -402,10 +402,12 @@ def test_box_mesh(ctx_getter, visualize=False): def test_as_python(): - from meshmode.mesh.generation import make_curve_mesh, cloverleaf - mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, 100), order=3) + from meshmode.mesh.generation import generate_box_mesh + mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),)) + # These implicitly compute these adjacency structures. mesh.nodal_adjacency + mesh.facial_adjacency_groups from meshmode.mesh import as_python code = as_python(mesh)