From 13b9cc32fcb023d6c063a7e331f7188b435a72a3 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 5 Nov 2018 17:50:15 -0600 Subject: [PATCH 01/12] made changes to mesh.io.GmshMeshReceiver to read in and mark physical boundary tags --- meshmode/mesh/io.py | 115 +++++++++++++++++++++++++++++++++++++++----- 1 file changed, 103 insertions(+), 12 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index d19f041d..ac535819 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -59,6 +59,11 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.element_types = None self.element_markers = None self.tags = None + self.tag_to_my_index = None # map name of tag to index in + # attr:: self.tags + self.gmsh_tag_index_to_mine = None + self.node_2_elts = None # list of lists, each node maps to + # a list of elements it is a member of if mesh_construction_kwargs is None: mesh_construction_kwargs = {} @@ -70,9 +75,11 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # Preallocation not done for performance, but to assign values at indices # in random order. self.points = [None] * count + self.node_2_elts = [None] * count def add_node(self, node_nr, point): self.points[node_nr] = point + self.node_2_elts[node_nr] = [] def finalize_nodes(self): self.points = np.array(self.points, dtype=np.float64) @@ -83,20 +90,45 @@ 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 + + # only physical tags are supported + if tag_numbers and len(tag_numbers) > 1: + tag_numbers = [tag_numbers[0]] self.element_markers[element_nr] = tag_numbers + # record this element in node to element map + for node in lexicographic_nodes: + self.node_2_elts[node].append(element_nr) + 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.tag_to_my_index is None: + self.tag_to_my_index = {} + if self.gmsh_tag_index_to_mine is None: + self.gmsh_tag_index_to_mine = {} + + if name not in self.tag_to_my_index: + # add tag if new + self.tag_to_my_index[name] = len(self.tags) + self.gmsh_tag_index_to_mine[index] = len(self.tags) + self.tags.append((name, index, dimension)) + else: + # ensure not trying to add different tags with same name + ndx = self.tag_to_my_index[name] + _, prev_index, prev_dim = self.tags[ndx] + if index != prev_index or dimension != prev_dim: + raise ValueError("Distinct tags with the same name") def finalize_tags(self): pass @@ -118,22 +150,22 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # {{{ build vertex numbering - vertex_index_gmsh_to_mine = {} + vertex_gmsh_index_to_mine = {} for el_vertices, el_type in 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) # }}} # {{{ 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 @@ -144,6 +176,37 @@ class GmshMeshReceiver(GmshMeshReceiverBase): bulk_el_types = set() + # get elements of dimension mesh_bulk_dim which contain elt + def get_super_elts(elt): + super_elts = [] + for pos_super in self.node_2_elts[self.element_nodes[elt][0]]: + if self.element_types[pos_super].dimensions != mesh_bulk_dim: + continue + + contains_each_node = True + for node in self.element_nodes[elt][1:]: + if pos_super in self.node_2_elts[node]: + contains_each_node = False + break + if contains_each_node: + super_elts.append(pos_super) + + return super_elts + + # populate tags from elements of small dimension to elements of + # full dimension (mesh_bulk_dim) + for e in range(len(self.element_types)): + if self.element_types[e].dimensions == mesh_bulk_dim: + continue + + super_elts = get_super_elts(e) + for se in super_elts: + for tag in self.element_markers[e]: + if tag not in self.element_markers[se]: + self.element_markers[se].append(tag) + + + element_index_mine_to_gmsh = {} for group_el_type, ngroup_elements in six.iteritems(el_type_hist): if group_el_type.dimensions != mesh_bulk_dim: continue @@ -158,14 +221,16 @@ 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] + n_elements = len(element_index_mine_to_gmsh) + element_index_mine_to_gmsh[n_elements] = element i += 1 @@ -215,11 +280,37 @@ class GmshMeshReceiver(GmshMeshReceiverBase): else: is_conforming = mesh_bulk_dim < 3 - return Mesh( + # prepare bdy tags for Mesh + bdy_tags = None + if self.tags is not None: + bdy_tags = [t[0] for t in self.tags] + + m = Mesh( vertices, groups, is_conforming=is_conforming, + boundary_tags = bdy_tags, **self.mesh_construction_kwargs) + # now mark boundaries with boundary tags + if len(self.tags) > 0: + for igrp in range(len(m.groups)): + bdry_fagrp = m.facial_adjacency_groups[igrp].get(None, None) + + if bdry_fagrp is None: + continue + + for ndx, e in enumerate(bdry_fagrp.elements): + gmsh_e = element_index_mine_to_gmsh[e] + for tag_number in self.element_markers[gmsh_e]: + name, _, _ = self.tags[self.gmsh_tag_index_to_mine[tag_number]] + btag_bit = m.boundary_tag_bit(name) + + neg = bdry_fagrp.neighbors[ndx] + assert(neg < 0) + bdry_fagrp.neighbors[ndx] = -((-neg) | btag_bit) + + return m + # }}} -- GitLab From 7a197da0fdefd7a88a0ce07fa26457390f56f5b4 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 5 Nov 2018 18:49:51 -0600 Subject: [PATCH 02/12] fixed flake8 errors --- meshmode/mesh/io.py | 22 ++++++++++------------ 1 file changed, 10 insertions(+), 12 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index ac535819..27c47848 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -59,11 +59,10 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.element_types = None self.element_markers = None self.tags = None - self.tag_to_my_index = None # map name of tag to index in - # attr:: self.tags + self.tag_to_my_index = None # map name of tag to index in attr:: self.tags self.gmsh_tag_index_to_mine = None - self.node_2_elts = None # list of lists, each node maps to - # a list of elements it is a member of + # list of lists, each node maps to a list of elements it is a member of + self.node_2_elts = None if mesh_construction_kwargs is None: mesh_construction_kwargs = {} @@ -102,7 +101,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase): tag_numbers = [tag_numbers[0]] self.element_markers[element_nr] = tag_numbers - # record this element in node to element map + # record this element in node to element map for node in lexicographic_nodes: self.node_2_elts[node].append(element_nr) @@ -122,8 +121,8 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # add tag if new self.tag_to_my_index[name] = len(self.tags) self.gmsh_tag_index_to_mine[index] = len(self.tags) - self.tags.append((name, index, dimension)) - else: + self.tags.append((name, index, dimension)) + else: # ensure not trying to add different tags with same name ndx = self.tag_to_my_index[name] _, prev_index, prev_dim = self.tags[ndx] @@ -205,7 +204,6 @@ class GmshMeshReceiver(GmshMeshReceiverBase): if tag not in self.element_markers[se]: self.element_markers[se].append(tag) - element_index_mine_to_gmsh = {} for group_el_type, ngroup_elements in six.iteritems(el_type_hist): if group_el_type.dimensions != mesh_bulk_dim: @@ -222,7 +220,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase): i = 0 for element, (el_vertices, el_nodes, el_type) in enumerate(zip( - self.element_vertices, self.element_nodes, self.element_types)): + self.element_vertices, self.element_nodes, self.element_types)): if el_type is not group_el_type: continue @@ -288,21 +286,21 @@ class GmshMeshReceiver(GmshMeshReceiverBase): m = Mesh( vertices, groups, is_conforming=is_conforming, - boundary_tags = bdy_tags, + boundary_tags=bdy_tags, **self.mesh_construction_kwargs) # now mark boundaries with boundary tags if len(self.tags) > 0: for igrp in range(len(m.groups)): bdry_fagrp = m.facial_adjacency_groups[igrp].get(None, None) - + if bdry_fagrp is None: continue for ndx, e in enumerate(bdry_fagrp.elements): gmsh_e = element_index_mine_to_gmsh[e] for tag_number in self.element_markers[gmsh_e]: - name, _, _ = self.tags[self.gmsh_tag_index_to_mine[tag_number]] + name = self.tags[self.gmsh_tag_index_to_mine[tag_number]][0] btag_bit = m.boundary_tag_bit(name) neg = bdry_fagrp.neighbors[ndx] -- GitLab From aa62b23f0efa0093e3e4f6c6008862328f192482 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 12 Nov 2018 12:37:11 -0600 Subject: [PATCH 03/12] Ensured tags exist before taking length --- meshmode/mesh/io.py | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 27c47848..c1352789 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -58,9 +58,9 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.elements = None self.element_types = None self.element_markers = None - self.tags = None - self.tag_to_my_index = None # map name of tag to index in attr:: self.tags - self.gmsh_tag_index_to_mine = None + self.tags = [] + self.tag_to_my_index = {} # map name of tag to index in attr:: self.tags + self.gmsh_tag_index_to_mine = {} # list of lists, each node maps to a list of elements it is a member of self.node_2_elts = None @@ -110,15 +110,8 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # May raise ValueError if try to add different tags with the same name def add_tag(self, name, index, dimension): - if self.tags is None: - self.tags = [] - if self.tag_to_my_index is None: - self.tag_to_my_index = {} - if self.gmsh_tag_index_to_mine is None: - self.gmsh_tag_index_to_mine = {} - + # add tag if new if name not in self.tag_to_my_index: - # add tag if new self.tag_to_my_index[name] = len(self.tags) self.gmsh_tag_index_to_mine[index] = len(self.tags) self.tags.append((name, index, dimension)) -- GitLab From ce345bc34d90cb159ea63b4647238dd277458e8d Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 12 Nov 2018 18:47:58 -0600 Subject: [PATCH 04/12] changed some names, made populating boundary tags to higher dimensional elements faster --- meshmode/mesh/io.py | 72 +++++++++++++++++++++++++++------------------ 1 file changed, 43 insertions(+), 29 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index c1352789..24bfa694 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -61,8 +61,8 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.tags = [] self.tag_to_my_index = {} # map name of tag to index in attr:: self.tags self.gmsh_tag_index_to_mine = {} - # list of lists, each node maps to a list of elements it is a member of - self.node_2_elts = None + # list of sets, each node maps to a set of elements it is a member of + self.node_to_containing_elements = None if mesh_construction_kwargs is None: mesh_construction_kwargs = {} @@ -74,11 +74,11 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # Preallocation not done for performance, but to assign values at indices # in random order. self.points = [None] * count - self.node_2_elts = [None] * count + self.node_to_containing_elements = [None] * count def add_node(self, node_nr, point): self.points[node_nr] = point - self.node_2_elts[node_nr] = [] + self.node_to_containing_elements[node_nr] = set() def finalize_nodes(self): self.points = np.array(self.points, dtype=np.float64) @@ -103,7 +103,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # record this element in node to element map for node in lexicographic_nodes: - self.node_2_elts[node].append(element_nr) + self.node_to_containing_elements[node].add(element_nr) def finalize_elements(self): pass @@ -168,34 +168,48 @@ class GmshMeshReceiver(GmshMeshReceiverBase): bulk_el_types = set() - # get elements of dimension mesh_bulk_dim which contain elt - def get_super_elts(elt): - super_elts = [] - for pos_super in self.node_2_elts[self.element_nodes[elt][0]]: - if self.element_types[pos_super].dimensions != mesh_bulk_dim: - continue - - contains_each_node = True - for node in self.element_nodes[elt][1:]: - if pos_super in self.node_2_elts[node]: - contains_each_node = False - break - if contains_each_node: - super_elts.append(pos_super) - - return super_elts + """Returns a set of the indices of elements with dimension mesh_bulk_dim + which contain all the nodes of elt. + + :arg element_ndx: The index of an element + """ + def get_higher_dim_element(element_ndx): + # Take the intersection of the sets of elements which + # contain at least one of the nodes used by the + # element at element_ndx + higher_dim_elts = None + for node in self.element_nodes[element_ndx]: + if higher_dim_elts is None: + higher_dim_elts = set(self.node_to_containing_elements[node]) + else: + higher_dim_elts &= self.node_to_containing_elements[node] + + # only keep elements of dimension mesh_bulk_dim + higher_dim_elts = {e for e in higher_dim_elts + if self.element_types[e].dimensions + == mesh_bulk_dim} + + # if no higher dimension elements, return empty set + if higher_dim_elts is None: + higher_dim_elts = set() + + return higher_dim_elts # populate tags from elements of small dimension to elements of # full dimension (mesh_bulk_dim) - for e in range(len(self.element_types)): - if self.element_types[e].dimensions == mesh_bulk_dim: - continue + if self.tags: + for elt_ndx in range(len(self.element_types)): + if self.element_types[elt_ndx].dimensions == mesh_bulk_dim: + continue + # if this element has no tags, continue + if not self.element_markers[elt_ndx]: + continue - super_elts = get_super_elts(e) - for se in super_elts: - for tag in self.element_markers[e]: - if tag not in self.element_markers[se]: - self.element_markers[se].append(tag) + 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) element_index_mine_to_gmsh = {} for group_el_type, ngroup_elements in six.iteritems(el_type_hist): -- GitLab From 66d6385232cc0c6667af00431e4ca70366d73364 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 12 Nov 2018 18:55:40 -0600 Subject: [PATCH 05/12] moved docstring to inside function --- meshmode/mesh/io.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 24bfa694..5f9f18c8 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -168,12 +168,12 @@ class GmshMeshReceiver(GmshMeshReceiverBase): bulk_el_types = set() - """Returns a set of the indices of elements with dimension mesh_bulk_dim - which contain all the nodes of elt. - - :arg element_ndx: The index of an element - """ def get_higher_dim_element(element_ndx): + """Returns a set of the indices of elements with dimension + mesh_bulk_dim which contain all the nodes of elt. + + :arg element_ndx: The index of an element + """ # Take the intersection of the sets of elements which # contain at least one of the nodes used by the # element at element_ndx -- GitLab From 784e06849ef23b58394aaded276da486f209b10a Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 12 Nov 2018 20:30:33 -0600 Subject: [PATCH 06/12] added test for boundary tags along with accompanying msh file --- meshmode/mesh/io.py | 4 +- test/annulus.msh | 210 ++++++++++++++++++++++++++++++++++++++++++ test/test_meshmode.py | 43 +++++++++ 3 files changed, 255 insertions(+), 2 deletions(-) create mode 100644 test/annulus.msh diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 5f9f18c8..750348f8 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -169,7 +169,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase): bulk_el_types = set() def get_higher_dim_element(element_ndx): - """Returns a set of the indices of elements with dimension + """Returns a set of the indices of elements with dimension mesh_bulk_dim which contain all the nodes of elt. :arg element_ndx: The index of an element @@ -297,7 +297,7 @@ class GmshMeshReceiver(GmshMeshReceiverBase): **self.mesh_construction_kwargs) # now mark boundaries with boundary tags - if len(self.tags) > 0: + if self.tags: for igrp in range(len(m.groups)): bdry_fagrp = m.facial_adjacency_groups[igrp].get(None, None) diff --git a/test/annulus.msh b/test/annulus.msh new file mode 100644 index 00000000..b105aa3f --- /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 e7ccf52a..ae99c829 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -76,6 +76,49 @@ 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 nbrs in 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) + +# }}} + + # {{{ convergence of boundary interpolation @pytest.mark.parametrize("group_factory", [ -- GitLab From 5351d2366b2ddf75a33a63f48bc365879423970b Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Tue, 13 Nov 2018 21:37:16 -0600 Subject: [PATCH 07/12] Modified Mesh class so that can pass in boundary tags during construction --- meshmode/mesh/__init__.py | 36 +++++++++++++++++++++++++- meshmode/mesh/io.py | 53 ++++++++++++++++++++------------------- 2 files changed, 62 insertions(+), 27 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 58f765c0..05d99f53 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -651,6 +651,7 @@ class Mesh(Record): nodal_adjacency=None, facial_adjacency_groups=None, boundary_tags=None, + element_boundary_tags=None, vertex_id_dtype=np.int32, element_id_dtype=np.int32, is_conforming=None): @@ -683,6 +684,25 @@ class Mesh(Record): :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. + :arg boundary_tags: One of two options: + *None*, in which case the mesh uses only default boundary + 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. + 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']] """ el_nr = 0 @@ -747,6 +767,10 @@ class Mesh(Record): element_id_dtype=np.dtype(element_id_dtype), is_conforming=is_conforming, ) + if is_conforming and element_boundary_tags and \ + facial_adjacency_groups is None: + self._facial_adjacency_groups = \ + _compute_facial_adjacency_from_vertices(self, element_boundary_tags) if not skip_tests: if node_vertex_consistency_tolerance is not False: @@ -996,7 +1020,11 @@ def _compute_nodal_adjacency_from_vertices(mesh): # {{{ vertex-based facial adjacency -def _compute_facial_adjacency_from_vertices(mesh): +def _compute_facial_adjacency_from_vertices(mesh, element_boundary_tags=None): + """ + :arg element_boundary_tags: See description of :arg:`element_boundary + tags1 in initializer of :class:`Mesh`. + """ # FIXME Native code would make this faster # create face_map, which is a mapping of @@ -1052,6 +1080,12 @@ def _compute_facial_adjacency_from_vertices(mesh): 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, diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 750348f8..e099d057 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -211,10 +211,23 @@ class GmshMeshReceiver(GmshMeshReceiverBase): 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 = {} 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 bulk_el_types.add(group_el_type) @@ -236,9 +249,21 @@ 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 @@ -285,37 +310,13 @@ class GmshMeshReceiver(GmshMeshReceiverBase): else: is_conforming = mesh_bulk_dim < 3 - # prepare bdy tags for Mesh - bdy_tags = None - if self.tags is not None: - bdy_tags = [t[0] for t in self.tags] - - m = Mesh( + return Mesh( vertices, groups, is_conforming=is_conforming, boundary_tags=bdy_tags, + element_boundary_tags=element_bdy_markers, **self.mesh_construction_kwargs) - # now mark boundaries with boundary tags - if self.tags: - for igrp in range(len(m.groups)): - bdry_fagrp = m.facial_adjacency_groups[igrp].get(None, None) - - if bdry_fagrp is None: - continue - - for ndx, e in enumerate(bdry_fagrp.elements): - gmsh_e = element_index_mine_to_gmsh[e] - for tag_number in self.element_markers[gmsh_e]: - name = self.tags[self.gmsh_tag_index_to_mine[tag_number]][0] - btag_bit = m.boundary_tag_bit(name) - - neg = bdry_fagrp.neighbors[ndx] - assert(neg < 0) - bdry_fagrp.neighbors[ndx] = -((-neg) | btag_bit) - - return m - # }}} -- GitLab From a1286ed28aca6cb90bdd294aaec5f63e33f00cbc Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Tue, 13 Nov 2018 21:37:29 -0600 Subject: [PATCH 08/12] fixed print formatting error --- test/test_meshmode.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index ae99c829..2d31af55 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -110,11 +110,11 @@ def test_boundary_tags(): # 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) + (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) + (num_marked_outer_bdy, num_on_outer_bdy)) # }}} -- GitLab From 5f7b24deaa4f4e4e598b36ab70f1dd68064f9923 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Tue, 13 Nov 2018 21:58:49 -0600 Subject: [PATCH 09/12] set facial_adjacency_groups to only compute once needed, even with boundary tags provided at construction of mesh --- meshmode/mesh/__init__.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 05d99f53..d9642e0d 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -717,6 +717,8 @@ class Mesh(Record): # {{{ boundary tags + self._element_boundary_tags = element_boundary_tags + if boundary_tags is None: boundary_tags = [] else: @@ -767,10 +769,6 @@ class Mesh(Record): element_id_dtype=np.dtype(element_id_dtype), is_conforming=is_conforming, ) - if is_conforming and element_boundary_tags and \ - facial_adjacency_groups is None: - self._facial_adjacency_groups = \ - _compute_facial_adjacency_from_vertices(self, element_boundary_tags) if not skip_tests: if node_vertex_consistency_tolerance is not False: @@ -888,7 +886,9 @@ class Mesh(Record): "be computed for known-conforming meshes") self._facial_adjacency_groups = \ - _compute_facial_adjacency_from_vertices(self) + _compute_facial_adjacency_from_vertices(self, + self._element_boundary_tags) + self._element_boundary_tags = None return self._facial_adjacency_groups -- GitLab From 40c7cf1ba9a44365c7118b3c0e3e34039f4ac2fa Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Mon, 26 Nov 2018 20:19:12 -0600 Subject: [PATCH 10/12] Fixed gmsh reading so that tagged boundary to appropriate face --- meshmode/mesh/__init__.py | 31 +++++------ meshmode/mesh/io.py | 112 ++++++++++++++++++++++++-------------- test/test_meshmode.py | 6 +- 3 files changed, 89 insertions(+), 60 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index d9642e0d..5f75653c 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 e099d057..cd7065c3 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 2d31af55..a9f0adbc 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']) + # }}} -- GitLab From 3832b3395c56dcf54b5a22e23a27dd3cb382a399 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Tue, 27 Nov 2018 21:30:12 -0600 Subject: [PATCH 11/12] Moved boundary tagging process out of __init__.py file and into io.py --- meshmode/mesh/__init__.py | 33 +--- meshmode/mesh/io.py | 312 +++++++++++++++++++++++--------------- 2 files changed, 188 insertions(+), 157 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 5f75653c..58f765c0 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -651,7 +651,6 @@ class Mesh(Record): nodal_adjacency=None, facial_adjacency_groups=None, boundary_tags=None, - element_boundary_tags=None, vertex_id_dtype=np.int32, element_id_dtype=np.int32, is_conforming=None): @@ -684,18 +683,6 @@ class Mesh(Record): :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. - :arg boundary_tags: One of two options: - *None*, in which case the mesh uses only default boundary - 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 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 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 @@ -710,8 +697,6 @@ class Mesh(Record): # {{{ boundary tags - self._element_boundary_tags = element_boundary_tags - if boundary_tags is None: boundary_tags = [] else: @@ -879,9 +864,7 @@ class Mesh(Record): "be computed for known-conforming meshes") self._facial_adjacency_groups = \ - _compute_facial_adjacency_from_vertices(self, - self._element_boundary_tags) - self._element_boundary_tags = None + _compute_facial_adjacency_from_vertices(self) return self._facial_adjacency_groups @@ -1013,11 +996,7 @@ def _compute_nodal_adjacency_from_vertices(mesh): # {{{ vertex-based facial adjacency -def _compute_facial_adjacency_from_vertices(mesh, element_boundary_tags=None): - """ - :arg element_boundary_tags: See description of :arg:`element_boundary - tags1 in initializer of :class:`Mesh`. - """ +def _compute_facial_adjacency_from_vertices(mesh): # FIXME Native code would make this faster # create face_map, which is a mapping of @@ -1136,14 +1115,6 @@ 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 cd7065c3..6696b0f4 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -58,11 +58,10 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.elements = None self.element_types = None self.element_markers = None - self.tags = [] - self.tag_to_my_index = {} # map name of tag to index in attr:: self.tags - self.gmsh_tag_index_to_mine = {} - # list of sets, each node maps to a set of elements it is a member of - self.node_to_containing_elements = None + self.tags = None + self.gmsh_tag_index_to_mine = None + # map set of face vertex indices to index of gmsh element + self.my_fvi_to_gmsh_elt_index = None if mesh_construction_kwargs is None: mesh_construction_kwargs = {} @@ -74,11 +73,9 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # Preallocation not done for performance, but to assign values at indices # in random order. self.points = [None] * count - self.node_to_containing_elements = [None] * count def add_node(self, node_nr, point): self.points[node_nr] = point - self.node_to_containing_elements[node_nr] = set() def finalize_nodes(self): self.points = np.array(self.points, dtype=np.float64) @@ -95,36 +92,183 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.element_vertices[element_nr] = vertex_nrs self.element_nodes[element_nr] = lexicographic_nodes self.element_types[element_nr] = element_type - - # only physical tags are supported if tag_numbers: - tag_numbers = [tag_numbers[0]] - self.element_markers[element_nr] = tag_numbers - - # record this element in node to element map - for node in lexicographic_nodes: - self.node_to_containing_elements[node].add(element_nr) + # 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): + 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 name not in self.tag_to_my_index: - self.tag_to_my_index[name] = len(self.tags) + if index not in self.gmsh_tag_index_to_mine: self.gmsh_tag_index_to_mine[index] = len(self.tags) - self.tags.append((name, index, dimension)) + self.tags.append((name, dimension)) else: - # ensure not trying to add different tags with same name - ndx = self.tag_to_my_index[name] - _, prev_index, prev_dim = self.tags[ndx] - if index != prev_index or dimension != prev_dim: - raise ValueError("Distinct tags with the same name") + # 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 + def _compute_facial_adjacency_groups(self, groups, boundary_tags, + element_id_dtype=np.int32, + face_id_dtype=np.int8): + if not groups: + return None + boundary_tag_to_index = {tag: i for i, tag in enumerate(boundary_tags)} + + def boundary_tag_bit(boundary_tag): + return 1 << boundary_tag_to_index[boundary_tag] + + # 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(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)) + + del igrp + del grp + + # 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") + + del face_tuples + 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(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=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) + element_faces.fill(-1) + neighbor_faces.fill(-1) + + neighbors.fill(-( + boundary_tag_bit(BTAG_ALL) + | boundary_tag_bit(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(groups)): + nb_count = group_count.get((igroup, ineighbor_group)) + if nb_count is not None: + 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) + element_faces.fill(-1) + neighbors.fill(-1) + neighbor_faces.fill(-1) + + grp_map[ineighbor_group] = FacialAdjacencyGroup( + igroup=igroup, + ineighbor_group=ineighbor_group, + elements=elements, + element_faces=element_faces, + neighbors=neighbors, + neighbor_faces=neighbor_faces) + + del igroup + del ineighbor_group + del grp_map + + # }}} + + # 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), (ineighbor_group, inb_el, inb_face) in [ + (face_tuples[0], face_tuples[1]), + (face_tuples[1], face_tuples[0]), + ]: + idx = fill_count.get((igroup, ineighbor_group), 0) + fill_count[igroup, ineighbor_group] = idx + 1 + + fagrp = facial_adjacency_groups[igroup][ineighbor_group] + 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((igroup, None), 0) + fill_count[igroup, None] = idx + 1 + + fagrp = facial_adjacency_groups[igroup][None] + fagrp.elements[idx] = iel + fagrp.element_faces[idx] = iface + # mark tags if present + if self.tags and self.my_fvi_to_gmsh_elt_index: + face_vertex_indices = groups[igroup].face_vertex_indices()[iface] + fvi = frozenset(groups[igroup].vertex_indices[ + iel, face_vertex_indices]) + gmsh_elt_index = self.my_fvi_to_gmsh_elt_index.get(fvi, None) + if gmsh_elt_index is not None: + tag = 0 + for t in self.element_markers[gmsh_elt_index]: + tag_name, _ = self.tags[self.gmsh_tag_index_to_mine[t]] + tag |= boundary_tag_bit(tag_name) + fagrp.neighbors[idx] = -(-(fagrp.neighbors[idx]) | tag) + + else: + raise RuntimeError("unexpected number of adjacent faces") + + return facial_adjacency_groups + + # }}} + def get_mesh(self): el_type_hist = {} for el_type in self.element_types: @@ -142,13 +286,16 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # {{{ build vertex numbering + self.my_fvi_to_gmsh_elt_index = {} vertex_gmsh_index_to_mine = {} - for el_vertices, el_type in zip( - self.element_vertices, self.element_types): + 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_gmsh_index_to_mine: vertex_gmsh_index_to_mine[gmsh_vertex_nr] = \ len(vertex_gmsh_index_to_mine) + el_grp_verts = {vertex_gmsh_index_to_mine[e] for e in el_vertices} + self.my_fvi_to_gmsh_elt_index[frozenset(el_grp_verts)] = element # }}} @@ -168,75 +315,9 @@ class GmshMeshReceiver(GmshMeshReceiverBase): bulk_el_types = set() - def get_higher_dim_element(element_ndx): - """Returns a set of the indices of elements with dimension - mesh_bulk_dim which contain all the nodes of elt. - - :arg element_ndx: The index of an element - """ - # Take the intersection of the sets of elements which - # contain at least one of the nodes used by the - # element at element_ndx - higher_dim_elts = None - for node in self.element_nodes[element_ndx]: - if higher_dim_elts is None: - higher_dim_elts = set(self.node_to_containing_elements[node]) - else: - higher_dim_elts &= self.node_to_containing_elements[node] - - # only keep elements of dimension mesh_bulk_dim - higher_dim_elts = {e for e in higher_dim_elts - if self.element_types[e].dimensions - == mesh_bulk_dim} - - # if no higher dimension elements, return empty set - if higher_dim_elts is None: - higher_dim_elts = set() - - return higher_dim_elts - - # 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)): - elt_markers = self.element_markers[elt_ndx] - if not elt_markers: - continue - 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 - # group-local element index to gmsh index - element_index_mine_to_gmsh = {} bulk_el_types.add(group_el_type) @@ -256,8 +337,6 @@ class GmshMeshReceiver(GmshMeshReceiverBase): nodes[:, i] = self.points[el_nodes].T vertex_indices[i] = [vertex_gmsh_index_to_mine[v_nr] for v_nr in el_vertices] - n_elements = len(element_index_mine_to_gmsh) - element_index_mine_to_gmsh[n_elements] = element i += 1 @@ -299,34 +378,6 @@ 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. @@ -335,16 +386,25 @@ class GmshMeshReceiver(GmshMeshReceiverBase): else: is_conforming = mesh_bulk_dim < 3 - # prepare bdy tags for Mesh - bdy_tags = None + # construct boundary tags for mesh + from meshmode.mesh import BTAG_ALL, BTAG_REALLY_ALL + boundary_tags = [BTAG_ALL, BTAG_REALLY_ALL] if self.tags: - bdy_tags = [tag for tag, _, dim in self.tags if dim == mesh_bulk_dim - 1] + 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: + facial_adjacency_groups = self._compute_facial_adjacency_groups( + groups, boundary_tags) return Mesh( vertices, groups, is_conforming=is_conforming, - boundary_tags=bdy_tags, - element_boundary_tags=group_elt_face_to_tags, + facial_adjacency_groups=facial_adjacency_groups, + boundary_tags=boundary_tags, **self.mesh_construction_kwargs) # }}} -- GitLab From 61b98f60d85535c7eec172c5a651575e9856df46 Mon Sep 17 00:00:00 2001 From: Ben Sepanski Date: Wed, 28 Nov 2018 11:14:00 -0600 Subject: [PATCH 12/12] removed duplicate code --- meshmode/mesh/__init__.py | 60 +++++++++---- meshmode/mesh/io.py | 174 ++++---------------------------------- 2 files changed, 62 insertions(+), 172 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 58f765c0..2c799478 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 6696b0f4..e4adf7ce 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -60,8 +60,6 @@ class GmshMeshReceiver(GmshMeshReceiverBase): self.element_markers = None self.tags = None self.gmsh_tag_index_to_mine = None - # map set of face vertex indices to index of gmsh element - self.my_fvi_to_gmsh_elt_index = None if mesh_construction_kwargs is None: mesh_construction_kwargs = {} @@ -120,155 +118,6 @@ class GmshMeshReceiver(GmshMeshReceiverBase): def finalize_tags(self): pass - def _compute_facial_adjacency_groups(self, groups, boundary_tags, - element_id_dtype=np.int32, - face_id_dtype=np.int8): - if not groups: - return None - boundary_tag_to_index = {tag: i for i, tag in enumerate(boundary_tags)} - - def boundary_tag_bit(boundary_tag): - return 1 << boundary_tag_to_index[boundary_tag] - - # 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(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)) - - del igrp - del grp - - # 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") - - del face_tuples - 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(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=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) - element_faces.fill(-1) - neighbor_faces.fill(-1) - - neighbors.fill(-( - boundary_tag_bit(BTAG_ALL) - | boundary_tag_bit(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(groups)): - nb_count = group_count.get((igroup, ineighbor_group)) - if nb_count is not None: - 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) - element_faces.fill(-1) - neighbors.fill(-1) - neighbor_faces.fill(-1) - - grp_map[ineighbor_group] = FacialAdjacencyGroup( - igroup=igroup, - ineighbor_group=ineighbor_group, - elements=elements, - element_faces=element_faces, - neighbors=neighbors, - neighbor_faces=neighbor_faces) - - del igroup - del ineighbor_group - del grp_map - - # }}} - - # 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), (ineighbor_group, inb_el, inb_face) in [ - (face_tuples[0], face_tuples[1]), - (face_tuples[1], face_tuples[0]), - ]: - idx = fill_count.get((igroup, ineighbor_group), 0) - fill_count[igroup, ineighbor_group] = idx + 1 - - fagrp = facial_adjacency_groups[igroup][ineighbor_group] - 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((igroup, None), 0) - fill_count[igroup, None] = idx + 1 - - fagrp = facial_adjacency_groups[igroup][None] - fagrp.elements[idx] = iel - fagrp.element_faces[idx] = iface - # mark tags if present - if self.tags and self.my_fvi_to_gmsh_elt_index: - face_vertex_indices = groups[igroup].face_vertex_indices()[iface] - fvi = frozenset(groups[igroup].vertex_indices[ - iel, face_vertex_indices]) - gmsh_elt_index = self.my_fvi_to_gmsh_elt_index.get(fvi, None) - if gmsh_elt_index is not None: - tag = 0 - for t in self.element_markers[gmsh_elt_index]: - tag_name, _ = self.tags[self.gmsh_tag_index_to_mine[t]] - tag |= boundary_tag_bit(tag_name) - fagrp.neighbors[idx] = -(-(fagrp.neighbors[idx]) | tag) - - else: - raise RuntimeError("unexpected number of adjacent faces") - - return facial_adjacency_groups - - # }}} - def get_mesh(self): el_type_hist = {} for el_type in self.element_types: @@ -286,7 +135,8 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # {{{ build vertex numbering - self.my_fvi_to_gmsh_elt_index = {} + # 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)): @@ -294,8 +144,17 @@ class GmshMeshReceiver(GmshMeshReceiverBase): 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) - el_grp_verts = {vertex_gmsh_index_to_mine[e] for e in el_vertices} - self.my_fvi_to_gmsh_elt_index[frozenset(el_grp_verts)] = element + 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 # }}} @@ -397,8 +256,10 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # compute facial adjacency for Mesh if there is tag information facial_adjacency_groups = None if is_conforming and self.tags: - facial_adjacency_groups = self._compute_facial_adjacency_groups( - groups, boundary_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, @@ -580,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, } -- GitLab