diff --git a/meshmode/interop/firedrake/mesh.py b/meshmode/interop/firedrake/mesh.py index bbc9e762f..48235dc6a 100644 --- a/meshmode/interop/firedrake/mesh.py +++ b/meshmode/interop/firedrake/mesh.py @@ -248,14 +248,14 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, # groups bdy_tags = _get_firedrake_boundary_tags( top, tag_induced_boundary=cells_to_use is not None) - boundary_tag_to_index = {bdy_tag: i for i, bdy_tag in enumerate(bdy_tags)} + from meshmode.mesh import make_tag_to_index, get_tag_bit + boundary_tag_to_index = make_tag_to_index(bdy_tags) marker_to_neighbor_value = {} - from meshmode.mesh import _boundary_tag_bit # for convenience, # None maps to the boundary tag for a boundary facet with no marker marker_to_neighbor_value[None] = \ - -(_boundary_tag_bit(bdy_tags, boundary_tag_to_index, BTAG_REALLY_ALL) - | _boundary_tag_bit(bdy_tags, boundary_tag_to_index, BTAG_ALL)) + -(get_tag_bit(boundary_tag_to_index, BTAG_REALLY_ALL) + | get_tag_bit(boundary_tag_to_index, BTAG_ALL)) # firedrake exterior facets with no marker are assigned the # a dummy marker from firedrake.mesh import unmarked as fd_unmarked @@ -263,7 +263,7 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, # Now figure out the appropriate tags for each firedrake markers for marker in top.exterior_facets.unique_markers: marker_to_neighbor_value[marker] = \ - -(_boundary_tag_bit(bdy_tags, boundary_tag_to_index, marker) + -(get_tag_bit(boundary_tag_to_index, marker) | -marker_to_neighbor_value[None]) # {{{ build the FacialAdjacencyGroup for internal connectivity @@ -314,12 +314,10 @@ def _get_firedrake_facial_adjacency_groups(fdrake_mesh_topology, newly_created_exterior_facs) new_ext_elements = int_elements[newly_created_exterior_facs] new_ext_element_faces = int_element_faces[newly_created_exterior_facs] - new_ext_neighbor_tag = -(_boundary_tag_bit(bdy_tags, - boundary_tag_to_index, - BTAG_REALLY_ALL) - | _boundary_tag_bit(bdy_tags, - boundary_tag_to_index, - BTAG_INDUCED_BOUNDARY)) + new_ext_neighbor_tag = -(get_tag_bit(boundary_tag_to_index, + BTAG_REALLY_ALL) + | get_tag_bit(boundary_tag_to_index, + BTAG_INDUCED_BOUNDARY)) new_ext_neighbors = np.full(new_ext_elements.shape, new_ext_neighbor_tag, dtype=IntType) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 583536863..0d2011a1e 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -38,10 +38,21 @@ .. autoclass:: FacialAdjacencyGroup .. autoclass:: InterPartitionAdjacencyGroup +.. autofunction:: make_region_tags +.. autofunction:: make_boundary_tags +.. autofunction:: make_tag_to_index +.. autofunction:: get_tag_bit + .. autofunction:: as_python .. autofunction:: check_bc_coverage .. autofunction:: is_boundary_tag_empty +Predefined Region tags +---------------------- + +.. autoclass:: RTAG_NONE +.. autoclass:: RTAG_ALL + Predefined Boundary tags ------------------------ @@ -56,12 +67,22 @@ # {{{ element tags +class RTAG_NONE: # noqa: N801 + """A region tag representing an empty region.""" + pass + + +class RTAG_ALL: # noqa: N801 + """A region tag representing all regions.""" + pass + + class BTAG_NONE: # noqa: N801 - """A boundary tag representing an empty boundary or volume.""" + """A boundary tag representing an empty boundary.""" class BTAG_ALL: # noqa: N801 - """A boundary tag representing the entire boundary or volume. + """A boundary tag representing the entire boundary. In the case of the boundary, :class:`BTAG_ALL` does not include rank boundaries, or, more generally, anything tagged with :class:`BTAG_NO_BOUNDARY`. @@ -134,7 +155,10 @@ class BTAG_INDUCED_BOUNDARY(BTAG_NO_BOUNDARY): # noqa: N801 # firedrakeproject.org seems to reject connections from Github. -SYSTEM_TAGS = {BTAG_NONE, BTAG_ALL, BTAG_REALLY_ALL, BTAG_NO_BOUNDARY, +SYSTEM_RTAGS = {RTAG_NONE, RTAG_ALL} + + +SYSTEM_BTAGS = {BTAG_NONE, BTAG_ALL, BTAG_REALLY_ALL, BTAG_NO_BOUNDARY, BTAG_PARTITION, BTAG_INDUCED_BOUNDARY} # }}} @@ -163,6 +187,11 @@ class MeshElementGroup(Record): *(dim, nunit_nodes)* + .. attribute:: regions + + An array *(nelements)* of integers, with the bits of ``regions[i]`` + indicating the mesh regions that contain element ``i``. + .. attribute:: element_nr_base Lowest element number in this element group. @@ -192,7 +221,7 @@ class MeshElementGroup(Record): """ def __init__(self, order, vertex_indices, nodes, - element_nr_base=None, node_nr_base=None, + regions=None, element_nr_base=None, node_nr_base=None, unit_nodes=None, dim=None, **kwargs): """ :arg order: the maximum total degree used for interpolation. @@ -210,6 +239,7 @@ def __init__(self, order, vertex_indices, nodes, vertex_indices=vertex_indices, nodes=nodes, unit_nodes=unit_nodes, + regions=regions, element_nr_base=element_nr_base, node_nr_base=node_nr_base, **kwargs) @@ -289,7 +319,7 @@ def __ne__(self, other): class _ModepyElementGroup(MeshElementGroup): def __init__(self, order, vertex_indices, nodes, - element_nr_base=None, node_nr_base=None, + regions=None, element_nr_base=None, node_nr_base=None, unit_nodes=None, dim=None, **kwargs): """ :arg order: the maximum total degree used for interpolation. @@ -341,6 +371,7 @@ def __init__(self, order, vertex_indices, nodes, f" got {vertex_indices.shape[-1]}") super().__init__(order, vertex_indices, nodes, + regions=regions, element_nr_base=element_nr_base, node_nr_base=node_nr_base, unit_nodes=unit_nodes, @@ -630,6 +661,15 @@ class Mesh(Record): (Note that element groups are not necessarily geometrically contiguous like the figure may suggest.) + .. attribute:: region_tags + + A list of region tag identifiers. :class:`RTAG_ALL` is guaranteed to exist. + + .. attribute:: rtag_to_index + + A mapping that maps region tag identifiers to their + corresponding index. + .. attribute:: boundary_tags A list of boundary tag identifiers. :class:`BTAG_ALL` and @@ -640,11 +680,6 @@ class Mesh(Record): A mapping that maps boundary tag identifiers to their corresponding index. - .. note:: - - Elements of :attr:`boundary_tags` that do not cover any - part of the boundary will not be keys in this dictionary. - .. attribute:: vertex_id_dtype .. attribute:: element_id_dtype @@ -667,6 +702,7 @@ def __init__(self, vertices, groups, *, skip_tests=False, skip_element_orientation_test=False, nodal_adjacency=None, facial_adjacency_groups=None, + region_tags=None, boundary_tags=None, vertex_id_dtype=np.int32, element_id_dtype=np.int32, @@ -710,20 +746,23 @@ def __init__(self, vertices, groups, *, skip_tests=False, el_nr += ng.nelements node_nr += ng.nnodes - # {{{ boundary tags + # {{{ region tags - if boundary_tags is None: - boundary_tags = [] - else: - boundary_tags = boundary_tags[:] + region_tags = make_region_tags(user_tags=region_tags) + + max_region_tag_count = int( + np.log(np.iinfo(element_id_dtype).max)/np.log(2)) + if len(region_tags) > max_region_tag_count: + raise ValueError("too few bits in element_id_dtype to represent all " + "region tags") + + rtag_to_index = make_tag_to_index(region_tags) + + # }}} + + # {{{ boundary tags - if BTAG_NONE in boundary_tags: - raise ValueError("BTAG_NONE is not allowed to be part of " - "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) + boundary_tags = make_boundary_tags(user_tags=boundary_tags) max_boundary_tag_count = int( np.log(np.iinfo(element_id_dtype).max)/np.log(2)) @@ -731,8 +770,7 @@ def __init__(self, vertices, groups, *, skip_tests=False, raise ValueError("too few bits in element_id_dtype to represent all " "boundary tags") - btag_to_index = { - btag: i for i, btag in enumerate(boundary_tags)} + btag_to_index = make_tag_to_index(boundary_tags) # }}} @@ -759,6 +797,8 @@ def __init__(self, vertices, groups, *, skip_tests=False, self, vertices=vertices, groups=new_groups, _nodal_adjacency=nodal_adjacency, _facial_adjacency_groups=facial_adjacency_groups, + region_tags=region_tags, + rtag_to_index=rtag_to_index, boundary_tags=boundary_tags, btag_to_index=btag_to_index, vertex_id_dtype=np.dtype(vertex_id_dtype), @@ -822,10 +862,10 @@ def set_if_not_present(name, from_name=None): set_if_not_present("vertices") if "groups" not in kwargs: kwargs["groups"] = [group.copy() for group in self.groups] + set_if_not_present("region_tags") set_if_not_present("boundary_tags") set_if_not_present("nodal_adjacency", "_nodal_adjacency") set_if_not_present("facial_adjacency_groups", "_facial_adjacency_groups") - set_if_not_present("boundary_tags") set_if_not_present("vertex_id_dtype") set_if_not_present("element_id_dtype") set_if_not_present("is_conforming") @@ -899,9 +939,11 @@ def facial_adjacency_groups(self): return self._facial_adjacency_groups + def region_tag_bit(self, region_tag): + return get_tag_bit(self.rtag_to_index, region_tag) + def boundary_tag_bit(self, boundary_tag): - return _boundary_tag_bit(self.boundary_tags, self.btag_to_index, - boundary_tag) + return get_tag_bit(self.btag_to_index, boundary_tag) def __eq__(self, other): return ( @@ -912,6 +954,7 @@ def __eq__(self, other): and self.element_id_dtype == other.element_id_dtype and self._nodal_adjacency == other._nodal_adjacency and self._facial_adjacency_groups == other._facial_adjacency_groups + and self.region_tags == other.region_tags and self.boundary_tags == other.boundary_tags and self.is_conforming == other.is_conforming) @@ -1029,17 +1072,55 @@ def _compute_nodal_adjacency_from_vertices(mesh): # }}} -# {{{ boundary tag to bit +# {{{ tags + +def make_region_tags(user_tags=None): + """Create a region tag list, optionally including extra *user_tags*.""" + region_tags = [] + + if user_tags is not None: + if RTAG_NONE in user_tags: + raise ValueError("RTAG_NONE is not allowed to be part of region_tags") + region_tags.extend(user_tags) + + if RTAG_ALL not in region_tags: + region_tags.append(RTAG_ALL) + + return region_tags + + +def make_boundary_tags(user_tags=None): + """Create a boundary tag list, optionally including extra *user_tags*.""" + boundary_tags = [] + + if user_tags is not None: + if BTAG_NONE in user_tags: + raise ValueError("BTAG_NONE is not allowed to be part of boundary_tags") + boundary_tags.extend(user_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) + + return boundary_tags + + +def make_tag_to_index(tags): + """Create a dict that maps tags to their respective index in the tag list.""" + return {tag: i for i, tag in enumerate(tags)} -def _boundary_tag_bit(boundary_tags, btag_to_index, boundary_tag): - if boundary_tag is BTAG_NONE: + +def get_tag_bit(tag_to_index, tag): + """Get the bit in a tag bitfield that corresponds to *tag*.""" + if tag is RTAG_NONE or tag is BTAG_NONE: return 0 - if boundary_tag not in boundary_tags: - raise ValueError("boundary tag '%s' is not known" % boundary_tag) + if tag not in tag_to_index: + raise ValueError("tag '%s' is not known" % tag) try: - return 1 << btag_to_index[boundary_tag] + return 1 << tag_to_index[tag] except KeyError: return 0 @@ -1135,10 +1216,10 @@ def _compute_facial_adjacency_from_vertices(groups, boundary_tags, if not groups: return None - boundary_tag_to_index = {tag: i for i, tag in enumerate(boundary_tags)} + boundary_tag_to_index = make_tag_to_index(boundary_tags) def boundary_tag_bit(boundary_tag): - return _boundary_tag_bit(boundary_tags, boundary_tag_to_index, boundary_tag) + return get_tag_bit(boundary_tag_to_index, boundary_tag) # Match up adjacent faces according to their vertex indices @@ -1411,13 +1492,27 @@ def check_bc_coverage(mesh, boundary_tags, incomplete_ok=False, # }}} +# {{{ is_region_tag_empty + +def is_region_tag_empty(mesh, region_tag): + """:returns: *True* if *region_tag* does not occur as part of *mesh*.""" + rtag_bit = mesh.region_tag_bit(region_tag) + if not rtag_bit: + return True + + for grp in mesh.groups: + if (grp.regions & rtag_bit).any(): + return False + + return True + +# }}} + + # {{{ is_boundary_tag_empty def is_boundary_tag_empty(mesh, boundary_tag): - """Return *True* if the corresponding boundary tag does not occur as part of - *mesh*. - """ - + """:returns: *True* if *boundary_tag* does not occur as part of *mesh*.""" btag_bit = mesh.boundary_tag_bit(boundary_tag) if not btag_bit: return True diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index aa3d78f75..b3595eaaf 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -311,12 +311,15 @@ def make_curve_mesh( t = t.ravel() nodes = curve_f(t).reshape(vertices.shape[0], nelements, -1) + regions = np.zeros(nelements, dtype=np.int32) + from meshmode.mesh import Mesh, SimplexElementGroup egroup = SimplexElementGroup( order, vertex_indices=vertex_indices, nodes=nodes, - unit_nodes=unit_nodes) + unit_nodes=unit_nodes, + regions=regions) mesh = Mesh( vertices=vertices, groups=[egroup], @@ -414,9 +417,11 @@ def make_group_from_vertices( # make contiguous nodes = nodes.copy() + regions = np.zeros(vertex_indices.shape[0], dtype=np.int32) + return group_cls( order, vertex_indices, nodes, - unit_nodes=unit_nodes) + unit_nodes=unit_nodes, regions=regions) # }}} diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 42ae557e3..34f024b05 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -51,7 +51,6 @@ def __init__(self, mesh_construction_kwargs=None): # Use data fields similar to meshpy.triangle.MeshInfo and # meshpy.tet.MeshInfo self.points = None - self.elements = None self.element_vertices = None self.element_nodes = None self.element_types = None @@ -141,8 +140,9 @@ def get_mesh(self): vertex_gmsh_index_to_mine[gmsh_vertex_nr] = \ len(vertex_gmsh_index_to_mine) if self.tags: + el_markers = self.element_markers[element] el_tag_indexes = [self.gmsh_tag_index_to_mine[t] for t in - self.element_markers[element]] + el_markers] if el_markers is not None else [] # record tags of boundary dimension el_tags = [self.tags[i][0] for i in el_tag_indexes if self.tags[i][1] == mesh_bulk_dim - 1] @@ -168,6 +168,21 @@ def get_mesh(self): from meshmode.mesh import (Mesh, SimplexElementGroup, TensorProductElementGroup) + # construct region tags for mesh + from meshmode.mesh import make_region_tags + gmsh_region_tags = [tag for tag, dim in self.tags if + dim == mesh_bulk_dim] if self.tags else None + region_tags = make_region_tags(user_tags=gmsh_region_tags) + + # construct boundary tags for mesh + from meshmode.mesh import make_boundary_tags + gmsh_boundary_tags = [tag for tag, dim in self.tags if + dim == mesh_bulk_dim-1] if self.tags else None + boundary_tags = make_boundary_tags(user_tags=gmsh_boundary_tags) + + from meshmode.mesh import make_tag_to_index, get_tag_bit + rtag_to_index = make_tag_to_index(region_tags) + bulk_el_types = set() for group_el_type, ngroup_elements in el_type_hist.items(): @@ -183,10 +198,11 @@ def get_mesh(self): vertex_indices = np.empty( (ngroup_elements, el_vertex_count), np.int32) + regions = np.zeros(ngroup_elements, np.int32) i = 0 - for el_vertices, el_nodes, el_type in zip( - self.element_vertices, self.element_nodes, self.element_types): + for element, (el_vertices, el_nodes, el_type) in enumerate(zip( + self.element_vertices, self.element_nodes, self.element_types)): if el_type is not group_el_type: continue @@ -195,6 +211,12 @@ def get_mesh(self): vertex_gmsh_index_to_mine[v_nr] for v_nr in el_vertices ] + el_markers = self.element_markers[element] + if el_markers is not None: + for t in el_markers: + regions[i] |= get_tag_bit(rtag_to_index, + self.tags[self.gmsh_tag_index_to_mine[t]][0]) + i += 1 import modepy as mp @@ -214,7 +236,8 @@ def get_mesh(self): group_el_type.order, vertex_indices, nodes, - unit_nodes=unit_nodes + unit_nodes=unit_nodes, + regions=regions ) if group.dim == 2: @@ -230,7 +253,8 @@ def get_mesh(self): group_el_type.order, vertex_indices[:, vertex_shuffle], nodes, - unit_nodes=unit_nodes + unit_nodes=unit_nodes, + regions=regions ) else: # NOTE: already checked above @@ -244,13 +268,6 @@ def get_mesh(self): else: is_conforming = mesh_bulk_dim < 3 - # construct boundary tags for mesh - from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL - boundary_tags = [BTAG_ALL, BTAG_REALLY_ALL] - if self.tags: - boundary_tags += [tag for tag, dim in self.tags if - dim == mesh_bulk_dim-1] - # compute facial adjacency for Mesh if there is tag information facial_adjacency_groups = None if is_conforming and self.tags: @@ -263,6 +280,7 @@ def get_mesh(self): vertices, groups, is_conforming=is_conforming, facial_adjacency_groups=facial_adjacency_groups, + region_tags=region_tags, boundary_tags=boundary_tags, **self.mesh_construction_kwargs) diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 695250a75..99536deb0 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -157,12 +157,19 @@ def _filter_mesh_groups(groups, selected_elements, vertex_id_dtype): # }}} - new_groups = [ - groups[i_old_grp].copy( - vertex_indices=new_vertex_indices[i_new_grp], - nodes=groups[i_old_grp].nodes[ - :, filtered_group_elements[i_new_grp], :].copy()) - for i_new_grp, i_old_grp in enumerate(new_group_to_old_group)] + new_groups = [] + for i_new_grp, i_old_grp in enumerate(new_group_to_old_group): + old_grp = groups[i_old_grp] + filtered_elems = filtered_group_elements[i_new_grp] + new_grp_vertex_indices = new_vertex_indices[i_new_grp] + new_grp_nodes = old_grp.nodes[:, filtered_elems, :].copy() + new_grp_regions = ( + old_grp.regions[filtered_elems].copy() + if old_grp.regions is not None else None) + new_groups.append(old_grp.copy( + vertex_indices=new_grp_vertex_indices, + nodes=new_grp_nodes, + regions=new_grp_regions)) return new_groups, group_to_new_group, required_vertex_indices @@ -587,17 +594,18 @@ def partition_mesh(mesh, part_per_element, part_num): group_neighbor_parts else set() boundary_tags = mesh.boundary_tags[:] - btag_to_index = {tag: i for i, tag in enumerate(boundary_tags)} - - def boundary_tag_bit(boundary_tag): - from meshmode.mesh import _boundary_tag_bit - return _boundary_tag_bit(boundary_tags, btag_to_index, boundary_tag) from meshmode.mesh import BTAG_PARTITION for i_neighbor_part in all_neighbor_parts: part_tag = BTAG_PARTITION(i_neighbor_part) boundary_tags.append(part_tag) - btag_to_index[part_tag] = len(boundary_tags)-1 + + from meshmode.mesh import make_tag_to_index + btag_to_index = make_tag_to_index(boundary_tags) + + def boundary_tag_bit(boundary_tag): + from meshmode.mesh import get_tag_bit + return get_tag_bit(btag_to_index, boundary_tag) inter_partition_adj_groups = _create_inter_partition_adjacency_groups(mesh, part_per_element, part_mesh_groups, all_neighbor_parts, @@ -614,6 +622,7 @@ def boundary_tag_bit(boundary_tag): part_mesh_groups, facial_adjacency_groups=part_facial_adj_groups, boundary_tags=boundary_tags, + region_tags=mesh.region_tags.copy(), is_conforming=mesh.is_conforming) return part_mesh, queried_elems