diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 58f765c0a821434e28b169456483e9d7d3bbf7f2..2c799478c31ded1f35b4c5f0c73cb26f6ca5ad2d 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -863,8 +863,11 @@ class Mesh(Record): raise DataUnavailable("facial_adjacency_groups can only " "be computed for known-conforming meshes") - self._facial_adjacency_groups = \ - _compute_facial_adjacency_from_vertices(self) + self._facial_adjacency_groups = _compute_facial_adjacency_from_vertices( + self.groups, + self.boundary_tags, + self.element_id_dtype, + self.face_id_dtype) return self._facial_adjacency_groups @@ -996,14 +999,27 @@ def _compute_nodal_adjacency_from_vertices(mesh): # {{{ vertex-based facial adjacency -def _compute_facial_adjacency_from_vertices(mesh): +def _compute_facial_adjacency_from_vertices(groups, boundary_tags, + element_id_dtype, + face_id_dtype, + face_vertex_indices_to_tags=None): + if not groups: + return None + boundary_tag_to_index = {tag: i for i, tag in enumerate(boundary_tags)} + + def boundary_tag_bit(boundary_tag): + try: + return 1 << boundary_tag_to_index[boundary_tag] + except KeyError: + raise 0 + # 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 igrp, grp in enumerate(groups): for fid, face_vertex_indices in enumerate(grp.face_vertex_indices()): all_fvi = grp.vertex_indices[:, face_vertex_indices] @@ -1031,18 +1047,19 @@ def _compute_facial_adjacency_from_vertices(mesh): del igrp # {{{ build facial_adjacency_groups data structure, still empty + from meshmode.mesh import FacialAdjacencyGroup, BTAG_ALL, BTAG_REALLY_ALL facial_adjacency_groups = [] - for igroup in range(len(mesh.groups)): + for igroup in range(len(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) + elements = np.empty(bdry_count, dtype=element_id_dtype) + element_faces = np.empty(bdry_count, dtype=face_id_dtype) + neighbors = np.empty(bdry_count, dtype=element_id_dtype) + neighbor_faces = np.zeros(bdry_count, dtype=face_id_dtype) # Ensure uninitialized entries get noticed elements.fill(-1) @@ -1050,8 +1067,8 @@ def _compute_facial_adjacency_from_vertices(mesh): neighbor_faces.fill(-1) neighbors.fill(-( - mesh.boundary_tag_bit(BTAG_ALL) - | mesh.boundary_tag_bit(BTAG_REALLY_ALL))) + boundary_tag_bit(BTAG_ALL) + | boundary_tag_bit(BTAG_REALLY_ALL))) grp_map[None] = FacialAdjacencyGroup( igroup=igroup, @@ -1061,13 +1078,13 @@ def _compute_facial_adjacency_from_vertices(mesh): neighbors=neighbors, neighbor_faces=neighbor_faces) - for ineighbor_group in range(len(mesh.groups)): + for ineighbor_group in range(len(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) + elements = np.empty(nb_count, dtype=element_id_dtype) + element_faces = np.empty(nb_count, dtype=face_id_dtype) + neighbors = np.empty(nb_count, dtype=element_id_dtype) + neighbor_faces = np.empty(nb_count, dtype=face_id_dtype) # Ensure uninitialized entries get noticed elements.fill(-1) @@ -1115,6 +1132,17 @@ def _compute_facial_adjacency_from_vertices(mesh): fagrp = facial_adjacency_groups[igroup][None] fagrp.elements[idx] = iel fagrp.element_faces[idx] = iface + # mark tags if present + if face_vertex_indices_to_tags: + face_vertex_indices = groups[igroup].face_vertex_indices()[iface] + fvi = frozenset(groups[igroup].vertex_indices[ + iel, face_vertex_indices]) + tags = face_vertex_indices_to_tags.get(fvi, None) + if tags is not None: + tag_mask = 0 + for tag in tags: + tag_mask |= boundary_tag_bit(tag) + fagrp.neighbors[idx] = -(-(fagrp.neighbors[idx]) | tag_mask) else: raise RuntimeError("unexpected number of adjacent faces") diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index d19f041df65de8f9bf6715dbc93437ae7dbbe919..e4adf7ce50983c863d82da66edfbc405c9c68f23 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -59,6 +59,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.element_types = None self.element_markers = None self.tags = None + self.gmsh_tag_index_to_mine = None if mesh_construction_kwargs is None: mesh_construction_kwargs = {} @@ -83,20 +84,36 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.element_nodes = [None] * count self.element_types = [None] * count self.element_markers = [None] * count - self.tags = [] def add_element(self, element_nr, element_type, vertex_nrs, lexicographic_nodes, tag_numbers): self.element_vertices[element_nr] = vertex_nrs self.element_nodes[element_nr] = lexicographic_nodes self.element_types[element_nr] = element_type - self.element_markers[element_nr] = tag_numbers + if tag_numbers: + # only physical tags are supported + physical_tag = tag_numbers[0] + self.element_markers[element_nr] = [physical_tag] def finalize_elements(self): pass + # May raise ValueError if try to add different tags with the same name def add_tag(self, name, index, dimension): - pass + if self.tags is None: + self.tags = [] + if self.gmsh_tag_index_to_mine is None: + self.gmsh_tag_index_to_mine = {} + # add tag if new + if index not in self.gmsh_tag_index_to_mine: + self.gmsh_tag_index_to_mine[index] = len(self.tags) + self.tags.append((name, dimension)) + else: + # ensure not trying to add different tags with same index + my_index = self.gmsh_tag_index_to_mine[index] + recorded_name, recorded_dim = self.tags[my_index] + if recorded_name != name or recorded_dim != dimension: + raise ValueError("Distinct tags with the same tag id") def finalize_tags(self): pass @@ -118,22 +135,35 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # {{{ build vertex numbering - vertex_index_gmsh_to_mine = {} - for el_vertices, el_type in zip( - self.element_vertices, self.element_types): + # map set of face vertex indices to list of tags associated to face + face_vertex_indices_to_tags = {} + vertex_gmsh_index_to_mine = {} + for element, (el_vertices, el_type) in enumerate(zip( + self.element_vertices, self.element_types)): for gmsh_vertex_nr in el_vertices: - if gmsh_vertex_nr not in vertex_index_gmsh_to_mine: - vertex_index_gmsh_to_mine[gmsh_vertex_nr] = \ - len(vertex_index_gmsh_to_mine) + if gmsh_vertex_nr not in vertex_gmsh_index_to_mine: + vertex_gmsh_index_to_mine[gmsh_vertex_nr] = \ + len(vertex_gmsh_index_to_mine) + if self.tags: + el_tag_indexes = [self.gmsh_tag_index_to_mine[t] for t in + self.element_markers[element]] + # 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] + el_grp_verts = {vertex_gmsh_index_to_mine[e] for e in el_vertices} + face_vertex_indices = frozenset(el_grp_verts) + if face_vertex_indices not in face_vertex_indices_to_tags: + face_vertex_indices_to_tags[face_vertex_indices] = [] + face_vertex_indices_to_tags[face_vertex_indices] += el_tags # }}} # {{{ build vertex array gmsh_vertex_indices, my_vertex_indices = \ - list(zip(*six.iteritems(vertex_index_gmsh_to_mine))) + list(zip(*six.iteritems(vertex_gmsh_index_to_mine))) vertices = np.empty( - (ambient_dim, len(vertex_index_gmsh_to_mine)), dtype=np.float64) + (ambient_dim, len(vertex_gmsh_index_to_mine)), dtype=np.float64) vertices[:, np.array(my_vertex_indices, np.intp)] = \ self.points[np.array(gmsh_vertex_indices, np.intp)].T @@ -158,13 +188,13 @@ class GmshMeshReceiver(GmshMeshReceiverBase): 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 nodes[:, i] = self.points[el_nodes].T - vertex_indices[i] = [vertex_index_gmsh_to_mine[v_nr] + vertex_indices[i] = [vertex_gmsh_index_to_mine[v_nr] for v_nr in el_vertices] i += 1 @@ -215,9 +245,27 @@ class GmshMeshReceiver(GmshMeshReceiverBase): 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] + boundary_tags = tuple(boundary_tags) + + # compute facial adjacency for Mesh if there is tag information + facial_adjacency_groups = None + if is_conforming and self.tags: + from meshmode.mesh import _compute_facial_adjacency_from_vertices + facial_adjacency_groups = _compute_facial_adjacency_from_vertices( + groups, boundary_tags, np.int32, np.int8, + face_vertex_indices_to_tags) + return Mesh( vertices, groups, is_conforming=is_conforming, + facial_adjacency_groups=facial_adjacency_groups, + boundary_tags=boundary_tags, **self.mesh_construction_kwargs) # }}} @@ -393,6 +441,7 @@ def to_json(mesh): "boundary_tags": [btag_to_json(btag) for btag in mesh.boundary_tags], "btag_to_index": dict( (btag_to_json(btag), value) + for btag, value in six.iteritems(mesh.btag_to_index)), "is_conforming": mesh.is_conforming, } diff --git a/test/annulus.msh b/test/annulus.msh new file mode 100644 index 0000000000000000000000000000000000000000..b105aa3fedfd1e801651f58ae6b4aa57c719a6d0 --- /dev/null +++ b/test/annulus.msh @@ -0,0 +1,210 @@ +$MeshFormat +2.2 0 8 +$EndMeshFormat +$PhysicalNames +3 +1 1 "outer_bdy" +1 2 "inner_bdy" +2 3 "surface" +$EndPhysicalNames +$Nodes +65 +1 0.5 0 0 +2 1 0 0 +3 0.442728012826605 0.2323615860218843 0 +4 0.284032373365578 0.4114919329468282 0 +5 0.06026834012766161 0.496354437049027 0 +6 -0.1773024435212677 0.4675081213427074 0 +7 -0.3742553740855504 0.3315613291203978 0 +8 -0.485470908713026 0.119657832143779 0 +9 -0.4854709087130261 -0.1196578321437787 0 +10 -0.3742553740855507 -0.3315613291203975 0 +11 -0.1773024435212679 -0.4675081213427074 0 +12 0.06026834012766116 -0.496354437049027 0 +13 0.2840323733655778 -0.4114919329468283 0 +14 0.4427280128266048 -0.2323615860218846 0 +15 0.970941817426052 0.2393156642875577 0 +16 0.8854560256532099 0.4647231720437685 0 +17 0.7485107481711012 0.6631226582407951 0 +18 0.5680647467311559 0.8229838658936564 0 +19 0.3546048870425358 0.9350162426854147 0 +20 0.1205366802553232 0.992708874098054 0 +21 -0.1205366802553229 0.992708874098054 0 +22 -0.3546048870425355 0.9350162426854148 0 +23 -0.5680647467311557 0.8229838658936565 0 +24 -0.7485107481711009 0.6631226582407955 0 +25 -0.8854560256532098 0.4647231720437687 0 +26 -0.9709418174260519 0.2393156642875581 0 +27 -1 1.224646799147353e-16 0 +28 -0.9709418174260521 -0.2393156642875574 0 +29 -0.8854560256532099 -0.4647231720437685 0 +30 -0.7485107481711013 -0.663122658240795 0 +31 -0.5680647467311559 -0.8229838658936564 0 +32 -0.3546048870425359 -0.9350162426854147 0 +33 -0.1205366802553236 -0.992708874098054 0 +34 0.1205366802553223 -0.9927088740980541 0 +35 0.3546048870425356 -0.9350162426854148 0 +36 0.5680647467311556 -0.8229838658936566 0 +37 0.7485107481711007 -0.6631226582407955 0 +38 0.8854560256532096 -0.4647231720437692 0 +39 0.9709418174260518 -0.2393156642875587 0 +40 0.6850831133807851 -0.1688578217853902 0 +41 -0.4008186261089604 0.5806860297800711 0 +42 -0.08642080595293153 -0.7117393709073221 0 +43 0.250204302333125 -0.659734524874673 0 +44 0.5366574057320725 -0.4754369744498491 0 +45 0.6961331124901092 0.1715814019523613 0 +46 0.5331079984486132 0.4722924739884488 0 +47 0.2542399547436574 0.6703756657375385 0 +48 -0.6348426322239812 0.3331911164962716 0 +49 -0.7169668665992207 2.289340507632594e-16 0 +50 -0.407283601489317 -0.5900521635914878 0 +51 -0.08642080595293102 0.7117393709073221 0 +52 -0.6348426322239811 -0.3331911164962713 0 +53 -0.7804257306297703 -0.1923576663407954 0 +54 -0.2850250891332481 -0.7515493938482497 0 +55 -0.601465483396546 0.5324187467308646 0 +56 0.8032623136056452 0.0003242357341630582 0 +57 0.7111027759033866 -0.3735817068465792 0 +58 0.4565718435150355 -0.6608877372876482 0 +59 0.4566297298396874 0.6617801944782985 0 +60 0.0965007009656307 -0.7974447092438742 0 +61 0.7119957035017039 0.3735315972877182 0 +62 0.09688523562135282 0.7979217028930111 0 +63 -0.7811982213823305 0.1925480682106987 0 +64 -0.2845375763037538 0.751178285934232 0 +65 -0.6030688547257194 -0.5342723842310627 0 +$EndNodes +$Elements +130 +1 1 2 2 1 1 3 +2 1 2 2 1 3 4 +3 1 2 2 1 4 5 +4 1 2 2 1 5 6 +5 1 2 2 1 6 7 +6 1 2 2 1 7 8 +7 1 2 2 1 8 9 +8 1 2 2 1 9 10 +9 1 2 2 1 10 11 +10 1 2 2 1 11 12 +11 1 2 2 1 12 13 +12 1 2 2 1 13 14 +13 1 2 2 1 14 1 +14 1 2 1 2 2 15 +15 1 2 1 2 15 16 +16 1 2 1 2 16 17 +17 1 2 1 2 17 18 +18 1 2 1 2 18 19 +19 1 2 1 2 19 20 +20 1 2 1 2 20 21 +21 1 2 1 2 21 22 +22 1 2 1 2 22 23 +23 1 2 1 2 23 24 +24 1 2 1 2 24 25 +25 1 2 1 2 25 26 +26 1 2 1 2 26 27 +27 1 2 1 2 27 28 +28 1 2 1 2 28 29 +29 1 2 1 2 29 30 +30 1 2 1 2 30 31 +31 1 2 1 2 31 32 +32 1 2 1 2 32 33 +33 1 2 1 2 33 34 +34 1 2 1 2 34 35 +35 1 2 1 2 35 36 +36 1 2 1 2 36 37 +37 1 2 1 2 37 38 +38 1 2 1 2 38 39 +39 1 2 1 2 39 2 +40 2 2 3 1 13 43 58 +41 2 2 3 1 7 55 48 +42 2 2 3 1 13 58 44 +43 2 2 3 1 5 47 62 +44 2 2 3 1 8 63 49 +45 2 2 3 1 4 59 47 +46 2 2 3 1 14 44 57 +47 2 2 3 1 12 42 60 +48 2 2 3 1 12 60 43 +49 2 2 3 1 3 61 46 +50 2 2 3 1 6 51 64 +51 2 2 3 1 10 52 65 +52 2 2 3 1 14 57 40 +53 2 2 3 1 3 45 61 +54 2 2 3 1 4 46 59 +55 2 2 3 1 7 41 55 +56 2 2 3 1 8 48 63 +57 2 2 3 1 10 65 50 +58 2 2 3 1 11 54 42 +59 2 2 3 1 11 50 54 +60 2 2 3 1 6 64 41 +61 2 2 3 1 9 49 53 +62 2 2 3 1 1 40 56 +63 2 2 3 1 5 62 51 +64 2 2 3 1 9 53 52 +65 2 2 3 1 1 56 45 +66 2 2 3 1 12 43 13 +67 2 2 3 1 1 14 40 +68 2 2 3 1 7 48 8 +69 2 2 3 1 4 47 5 +70 2 2 3 1 11 42 12 +71 2 2 3 1 3 46 4 +72 2 2 3 1 13 44 14 +73 2 2 3 1 8 49 9 +74 2 2 3 1 10 50 11 +75 2 2 3 1 1 45 3 +76 2 2 3 1 5 51 6 +77 2 2 3 1 9 52 10 +78 2 2 3 1 6 41 7 +79 2 2 3 1 35 43 60 +80 2 2 3 1 39 56 40 +81 2 2 3 1 33 42 54 +82 2 2 3 1 23 41 64 +83 2 2 3 1 27 49 63 +84 2 2 3 1 29 52 53 +85 2 2 3 1 15 61 45 +86 2 2 3 1 23 55 41 +87 2 2 3 1 27 53 49 +88 2 2 3 1 37 44 58 +89 2 2 3 1 33 60 42 +90 2 2 3 1 35 58 43 +91 2 2 3 1 37 57 44 +92 2 2 3 1 15 45 56 +93 2 2 3 1 17 46 61 +94 2 2 3 1 31 50 65 +95 2 2 3 1 17 59 46 +96 2 2 3 1 21 64 51 +97 2 2 3 1 25 48 55 +98 2 2 3 1 21 51 62 +99 2 2 3 1 29 65 52 +100 2 2 3 1 31 54 50 +101 2 2 3 1 19 47 59 +102 2 2 3 1 25 63 48 +103 2 2 3 1 19 62 47 +104 2 2 3 1 39 40 57 +105 2 2 3 1 2 56 39 +106 2 2 3 1 34 35 60 +107 2 2 3 1 27 28 53 +108 2 2 3 1 26 27 63 +109 2 2 3 1 24 25 55 +110 2 2 3 1 37 38 57 +111 2 2 3 1 30 31 65 +112 2 2 3 1 32 33 54 +113 2 2 3 1 15 16 61 +114 2 2 3 1 28 29 53 +115 2 2 3 1 16 17 61 +116 2 2 3 1 22 23 64 +117 2 2 3 1 33 34 60 +118 2 2 3 1 2 15 56 +119 2 2 3 1 23 24 55 +120 2 2 3 1 25 26 63 +121 2 2 3 1 36 37 58 +122 2 2 3 1 18 19 59 +123 2 2 3 1 19 20 62 +124 2 2 3 1 31 32 54 +125 2 2 3 1 17 18 59 +126 2 2 3 1 35 36 58 +127 2 2 3 1 38 39 57 +128 2 2 3 1 20 21 62 +129 2 2 3 1 21 22 64 +130 2 2 3 1 29 30 65 +$EndElements diff --git a/test/test_meshmode.py b/test/test_meshmode.py index e7ccf52afa6a321e2b5149ec55c725a05787930d..a9f0adbca854fabcc54ee5c4a71b44e8b544becb 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -76,6 +76,53 @@ def test_circle_mesh(do_plot=False): # }}} +# {{{ test boundary tags + +def test_boundary_tags(): + from meshmode.mesh.io import read_gmsh + # ensure tags are read in + mesh = read_gmsh('annulus.msh') + if not {'outer_bdy', 'inner_bdy'} <= set(mesh.boundary_tags): + print("Mesh boundary tags:", mesh.boundary_tags) + raise ValueError('Tags not saved by mesh') + + # correct answers + num_on_outer_bdy = 26 + num_on_inner_bdy = 13 + + # check how many elements are marked on each boundary + num_marked_outer_bdy = 0 + num_marked_inner_bdy = 0 + outer_btag_bit = mesh.boundary_tag_bit('outer_bdy') + inner_btag_bit = mesh.boundary_tag_bit('inner_bdy') + for igrp in range(len(mesh.groups)): + bdry_fagrp = mesh.facial_adjacency_groups[igrp].get(None, None) + + if bdry_fagrp is None: + continue + + for i, nbrs in enumerate(bdry_fagrp.neighbors): + if (-nbrs) & outer_btag_bit: + num_marked_outer_bdy += 1 + if (-nbrs) & inner_btag_bit: + num_marked_inner_bdy += 1 + + # raise errors if wrong number of elements marked + if num_marked_inner_bdy != num_on_inner_bdy: + raise ValueError("%i marked on inner boundary, should be %i" % + (num_marked_inner_bdy, num_on_inner_bdy)) + + if num_marked_outer_bdy != num_on_outer_bdy: + 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']) + +# }}} + + # {{{ convergence of boundary interpolation @pytest.mark.parametrize("group_factory", [