diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index d9642e0d2b63ea8e3e8345a3dc40ceaa26570784..5f75653c948f65b683a113132f84f522ae8bcb4e 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -689,20 +689,13 @@ class Mesh(Record): tags. A list of names of boundary tags, in which case the mesh appends any default tags not already included. :arg element_boundary_tags: - If facial_adjacency_groups is included or is_conforming is - False, nothing is done. + If *facial_adjacency_groups* is passed as an argument + or *is_conforming* is *False*, nothing is done. Else, one of two options: *None*, in which case the mesh marks only default boundary - tags. Else, a list of lists indicated which boundary tags - are associated to which elements. The [igrp] entry is a - list of the nonempty boundary tags associated to elements in group - at index igrp. Relative order of elements is preserved. - EXAMPLE: if group igrp has 3 elements, the first is tagged - as 'inner_bdy', the second has no boundary tags, and the - third is tagged 'outer_bdy' and 'exterior', then we would - have - element_boundary_tags[igrp] = - [['inner_bdy'], ['outer_bdy', 'exterior']] + tags. Else, a map which maps the index of the group, a + group-local index, and a face index to a list of tag names. + i.e. (igrp, igrp_elt, fid) to ['tag_one',...,'tag_n'] """ el_nr = 0 @@ -1080,12 +1073,6 @@ def _compute_facial_adjacency_from_vertices(mesh, element_boundary_tags=None): neighbors.fill(-( mesh.boundary_tag_bit(BTAG_ALL) | mesh.boundary_tag_bit(BTAG_REALLY_ALL))) - if element_boundary_tags: - for ndx, boundary_tags in enumerate(element_boundary_tags[igroup]): - tag = 0 - for bt in boundary_tags: - tag |= mesh.boundary_tag_bit(bt) - neighbors[ndx] = -((-neighbors[ndx]) | tag) grp_map[None] = FacialAdjacencyGroup( igroup=igroup, @@ -1149,6 +1136,14 @@ def _compute_facial_adjacency_from_vertices(mesh, element_boundary_tags=None): fagrp = facial_adjacency_groups[igroup][None] fagrp.elements[idx] = iel fagrp.element_faces[idx] = iface + # mark tags if present + if element_boundary_tags: + tags = element_boundary_tags.get(face_tuples[0], None) + if tags: + tag = 0 + for t in tags: + tag |= mesh.boundary_tag_bit(t) + fagrp.neighbors[idx] = -(-(fagrp.neighbors[idx]) | tag) else: raise RuntimeError("unexpected number of adjacent faces") diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index e099d0570a54509f2b5400ed17fba8a2859388f1..cd7065c36f9ec0801843ef64c99705d590cbdd5e 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -97,7 +97,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.element_types[element_nr] = element_type # only physical tags are supported - if tag_numbers and len(tag_numbers) > 1: + if tag_numbers: tag_numbers = [tag_numbers[0]] self.element_markers[element_nr] = tag_numbers @@ -195,39 +195,48 @@ class GmshMeshReceiver(GmshMeshReceiverBase): return higher_dim_elts - # populate tags from elements of small dimension to elements of - # full dimension (mesh_bulk_dim) + # elt_ndx maps to (dim of tag, face indices, tag) + elt_gmsh_index_to_tags = {} + # Record tags if self.tags: for elt_ndx in range(len(self.element_types)): - if self.element_types[elt_ndx].dimensions == mesh_bulk_dim: + elt_markers = self.element_markers[elt_ndx] + if not elt_markers: continue - # if this element has no tags, continue - if not self.element_markers[elt_ndx]: - continue - - higher_dim_elements = get_higher_dim_element(elt_ndx) - for higher_dim_elt in higher_dim_elements: - for tag in self.element_markers[elt_ndx]: - if tag not in self.element_markers[higher_dim_elt]: - self.element_markers[higher_dim_elt].append(tag) - - # prepare bdy tags for Mesh - bdy_tags = None - if self.tags: - bdy_tags = [t[0] for t in self.tags if t[-1] == mesh_bulk_dim - 1] - - # for each group, a list of non-empty boundary tags - element_bdy_markers = None - if self.tags: - element_bdy_markers = [] - - element_index_mine_to_gmsh = {} + elt_dim = self.element_types[elt_ndx].dimensions + + # if element is of dimension mesh_bulk_dim, record any tags + if elt_dim == mesh_bulk_dim: + face_vertices = tuple(range(len(self.element_vertices[elt_ndx]))) + if elt_ndx not in elt_gmsh_index_to_tags: + elt_gmsh_index_to_tags[elt_ndx] = [] + for tag_nr in elt_markers: + tag_name, _, dim = self.tags[ + self.gmsh_tag_index_to_mine[tag_nr]] + elt_gmsh_index_to_tags[elt_ndx].append( + (dim, face_vertices, tag_name)) + + # else, record the tag in the higher dimension elements + # with this as a face + else: + elt_vertices = set(self.element_vertices[elt_ndx]) + higher_dim_elements = get_higher_dim_element(elt_ndx) + for super_elt in higher_dim_elements: + # record tags + if super_elt not in elt_gmsh_index_to_tags: + elt_gmsh_index_to_tags[super_elt] = [] + for tag_nr in elt_markers: + tag_name, _, dim = self.tags[ + self.gmsh_tag_index_to_mine[tag_nr]] + elt_gmsh_index_to_tags[super_elt].append( + (dim, elt_vertices, tag_name)) + + group_elt_face_to_tags = {} for group_el_type, ngroup_elements in six.iteritems(el_type_hist): if group_el_type.dimensions != mesh_bulk_dim: continue - if self.tags: - boundary_markers = [] # list of nonempty boundary tags in - # preserving relative order of index of corresponding faces + # group-local element index to gmsh index + element_index_mine_to_gmsh = {} bulk_el_types.add(group_el_type) @@ -249,21 +258,9 @@ class GmshMeshReceiver(GmshMeshReceiverBase): for v_nr in el_vertices] n_elements = len(element_index_mine_to_gmsh) element_index_mine_to_gmsh[n_elements] = element - # record the tags associated to this element if any - if self.tags and self.element_markers[element]: - element_tags = [] - for t in self.element_markers[element]: - tag = self.tags[self.gmsh_tag_index_to_mine[t]] - if tag[-1] != mesh_bulk_dim - 1: - continue - element_tags.append(tag[0]) - if element_tags: - boundary_markers.append(element_tags) i += 1 - if element_bdy_markers is not None: - element_bdy_markers.append(boundary_markers) unit_nodes = (np.array(group_el_type.lexicographic_node_tuples(), dtype=np.float64).T/group_el_type.order)*2 - 1 @@ -302,6 +299,34 @@ class GmshMeshReceiver(GmshMeshReceiverBase): raise NotImplementedError("gmsh element type: %s" % type(group_el_type).__name__) + # record tags in group_elt_face_to_tags + if self.tags: + igrp = len(groups) + for elt_ndx in range(group.nelements): + gmsh_ndx = element_index_mine_to_gmsh[elt_ndx] + if gmsh_ndx not in elt_gmsh_index_to_tags: + continue + tags = elt_gmsh_index_to_tags[gmsh_ndx] + if tags: + for dim, elt_verts, name in tags: + # only record tags of boundary dimension + if dim != mesh_bulk_dim - 1: + continue + elt_verts = {vertex_gmsh_index_to_mine[e] + for e in elt_verts} + face_ndx = -1 + for fid, fvind in enumerate(group.face_vertex_indices()): + face_verts = {group.vertex_indices[elt_ndx][i] + for i in fvind} + if elt_verts == face_verts: + face_ndx = fid + break + assert face_ndx >= 0 + grp_el_face = (igrp, elt_ndx, face_ndx) + if grp_el_face not in group_elt_face_to_tags: + group_elt_face_to_tags[grp_el_face] = [] + group_elt_face_to_tags[grp_el_face].append(name) + groups.append(group) # FIXME: This is heuristic. @@ -310,11 +335,16 @@ class GmshMeshReceiver(GmshMeshReceiverBase): else: is_conforming = mesh_bulk_dim < 3 + # prepare bdy tags for Mesh + bdy_tags = None + if self.tags: + bdy_tags = [tag for tag, _, dim in self.tags if dim == mesh_bulk_dim - 1] + return Mesh( vertices, groups, is_conforming=is_conforming, boundary_tags=bdy_tags, - element_boundary_tags=element_bdy_markers, + element_boundary_tags=group_elt_face_to_tags, **self.mesh_construction_kwargs) # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 2d31af55d4aacfc3ad7293b6668cc34729814bb8..a9f0adbca854fabcc54ee5c4a71b44e8b544becb 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -101,7 +101,7 @@ def test_boundary_tags(): if bdry_fagrp is None: continue - for nbrs in bdry_fagrp.neighbors: + for i, nbrs in enumerate(bdry_fagrp.neighbors): if (-nbrs) & outer_btag_bit: num_marked_outer_bdy += 1 if (-nbrs) & inner_btag_bit: @@ -116,6 +116,10 @@ def test_boundary_tags(): raise ValueError("%i marked on outer boundary, should be %i" % (num_marked_outer_bdy, num_on_outer_bdy)) + # ensure boundary is covered + from meshmode.mesh import check_bc_coverage + check_bc_coverage(mesh, ['inner_bdy', 'outer_bdy']) + # }}}