From 3cd564d8ad1c6dba6c5fc0193f7241bc872555b7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 18 Apr 2018 11:54:45 -0500 Subject: [PATCH] Track whether meshes are known to be conforming --- meshmode/mesh/__init__.py | 38 ++++++++++++-- meshmode/mesh/generation.py | 21 +++----- meshmode/mesh/io.py | 31 ++++++++---- meshmode/mesh/processing.py | 24 ++++++--- meshmode/mesh/refinement/__init__.py | 41 ++++++++++----- meshmode/mesh/{ => refinement}/tesselate.py | 55 +++++++++++++++++++-- test/test_meshmode.py | 7 +-- test/test_refinement.py | 4 ++ 8 files changed, 168 insertions(+), 53 deletions(-) rename meshmode/mesh/{ => refinement}/tesselate.py (57%) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index d3368b44..89e97a13 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -523,6 +523,12 @@ class Mesh(Record): .. attribute:: element_id_dtype + .. attribute:: is_conforming + + *True* if it is known that all element interfaces are conforming. + *False* if it is known that some element interfaces are non-conforming. + *None* if neither of the two is known. + .. automethod:: __eq__ .. automethod:: __ne__ """ @@ -531,11 +537,12 @@ class Mesh(Record): def __init__(self, vertices, groups, skip_tests=False, node_vertex_consistency_tolerance=None, - nodal_adjacency=False, - facial_adjacency_groups=False, + nodal_adjacency=None, + facial_adjacency_groups=None, boundary_tags=None, vertex_id_dtype=np.int32, - element_id_dtype=np.int32): + element_id_dtype=np.int32, + is_conforming=None): """ The following are keyword-only: @@ -599,6 +606,12 @@ class Mesh(Record): # }}} + if is_conforming is not True: + if nodal_adjacency is None: + nodal_adjacency = False + if facial_adjacency_groups is None: + facial_adjacency_groups = False + if nodal_adjacency is not False and nodal_adjacency is not None: if not isinstance(nodal_adjacency, NodalAdjacency): nb_starts, nbs = nodal_adjacency @@ -617,6 +630,7 @@ class Mesh(Record): btag_to_index=btag_to_index, vertex_id_dtype=np.dtype(vertex_id_dtype), element_id_dtype=np.dtype(element_id_dtype), + is_conforming=is_conforming, ) if not skip_tests: @@ -704,7 +718,13 @@ class Mesh(Record): if self._nodal_adjacency is False: from meshmode import DataUnavailable raise DataUnavailable("nodal_adjacency") + elif self._nodal_adjacency is None: + if self.is_conforming is not True: + from meshmode import DataUnavailable + raise DataUnavailable("nodal_adjacency can only " + "be computed for known-conforming meshes") + self._nodal_adjacency = _compute_nodal_adjacency_from_vertices(self) return self._nodal_adjacency @@ -721,7 +741,13 @@ class Mesh(Record): if self._facial_adjacency_groups is False: from meshmode import DataUnavailable raise DataUnavailable("facial_adjacency_groups") + elif self._facial_adjacency_groups is None: + if self.is_conforming is not True: + from meshmode import DataUnavailable + raise DataUnavailable("facial_adjacency_groups can only " + "be computed for known-conforming meshes") + self._facial_adjacency_groups = \ _compute_facial_adjacency_from_vertices(self) @@ -744,7 +770,8 @@ class Mesh(Record): == other._nodal_adjacency) and (self._facial_adjacency_groups == other._facial_adjacency_groups) - and self.boundary_tags == other.boundary_tags) + and self.boundary_tags == other.boundary_tags + and self.is_conforming == other.is_conforming) def __ne__(self, other): return not self.__eq__(other) @@ -1082,7 +1109,8 @@ def as_python(mesh, function_name="make_mesh"): cg(" nodal_adjacency=%s," % el_con_str) cg(" facial_adjacency_groups=facial_adjacency_groups,") - cg(" boundary_tags=[%s])" % btags_str) + cg(" boundary_tags=[%s]," % btags_str) + cg(" is_conforming=%s)" % repr(mesh.is_conforming)) # FIXME: Handle facial adjacency, boundary tags diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 198c45f4..6f5d283c 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -281,9 +281,8 @@ def make_curve_mesh(curve_f, element_boundaries, order, mesh = Mesh( vertices=vertices, groups=[egroup], - nodal_adjacency=None, - facial_adjacency_groups=None, - node_vertex_consistency_tolerance=node_vertex_consistency_tolerance) + node_vertex_consistency_tolerance=node_vertex_consistency_tolerance, + is_conforming=True) if return_parametrization_points: return mesh, t @@ -420,8 +419,7 @@ def generate_icosahedron(r, order): from meshmode.mesh import Mesh return Mesh( vertices, [grp], - nodal_adjacency=None, - facial_adjacency_groups=None) + is_conforming=True) # }}} @@ -439,8 +437,7 @@ def generate_icosphere(r, order): from meshmode.mesh import Mesh return Mesh( mesh.vertices, [grp], - nodal_adjacency=None, - facial_adjacency_groups=None) + is_conforming=True) # }}} @@ -505,8 +502,7 @@ def generate_torus_and_cycle_vertices(r_outer, r_inner, return ( Mesh( vertices, [grp.copy(nodes=nodes)], - nodal_adjacency=None, - facial_adjacency_groups=None), + is_conforming=True), [idx(i, 0) for i in range(n_outer)], [idx(0, j) for j in range(n_inner)]) @@ -564,9 +560,9 @@ def refine_mesh_and_get_urchin_warper(order, m, n, est_rel_interp_tolerance, return Mesh( map_coords(mesh.vertices), groups, - nodal_adjacency=None, node_vertex_consistency_tolerance=False, - facial_adjacency_groups=None) + is_conforming=mesh.is_conforming, + ) unwarped_mesh = generate_icosphere(1, order=order) @@ -730,8 +726,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64, from meshmode.mesh import Mesh return Mesh(vertices, [grp], - nodal_adjacency=None, - facial_adjacency_groups=None) + is_conforming=True) # }}} diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index fd8a82a2..c663a3eb 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -137,10 +137,14 @@ class GmshMeshReceiver(GmshMeshReceiverBase): from meshmode.mesh import (Mesh, SimplexElementGroup, TensorProductElementGroup) + bulk_el_types = set() + for group_el_type, ngroup_elements in six.iteritems(el_type_hist): if group_el_type.dimensions != mesh_bulk_dim: continue + bulk_el_types.add(group_el_type) + nodes = np.empty((ambient_dim, ngroup_elements, el_type.node_count()), np.float64) el_vertex_count = group_el_type.vertex_count() @@ -198,15 +202,17 @@ class GmshMeshReceiver(GmshMeshReceiverBase): raise NotImplementedError("gmsh element type: %s" % type(group_el_type).__name__) - # Gmsh seems to produce elements in the opposite orientation - # of what we like. Flip them all. - groups.append(group) + # FIXME: This is heuristic. + if len(bulk_el_types) == 1: + is_conforming = True + else: + is_conforming = mesh_bulk_dim < 3 + return Mesh( vertices, groups, - nodal_adjacency=None, - facial_adjacency_groups=None) + is_conforming=is_conforming) # }}} @@ -287,8 +293,7 @@ def from_meshpy(mesh_info, order=1): return Mesh( vertices=vertices, groups=[grp], - nodal_adjacency=None, - facial_adjacency_groups=None) + is_conforming=True) # }}} @@ -320,8 +325,7 @@ def from_vertices_and_simplices(vertices, simplices, order=1, fix_orientation=Fa return Mesh( vertices=vertices, groups=[grp], - nodal_adjacency=None, - facial_adjacency_groups=None) + is_conforming=True) # }}} @@ -364,7 +368,13 @@ def to_json(mesh): } return { - "version": 0, + # VERSION 0: + # - initial version + # + # VERSION 1: + # - added is_conforming + + "version": 1, "vertices": mesh.vertices.tolist(), "groups": [group_to_json(group) for group in mesh.groups], "nodal_adjacency": nodal_adjacency_to_json(mesh), @@ -374,6 +384,7 @@ def to_json(mesh): "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/meshmode/mesh/processing.py b/meshmode/mesh/processing.py index 664355e0..36c38eda 100644 --- a/meshmode/mesh/processing.py +++ b/meshmode/mesh/processing.py @@ -137,11 +137,16 @@ def partition_mesh(mesh, part_per_element, part_nr): if group_nr not in skip_groups: mesh_group = mesh.groups[group_nr] new_mesh_groups.append( - type(mesh_group)(mesh_group.order, new_indices[group_nr], - new_nodes[group_nr], unit_nodes=mesh_group.unit_nodes)) + type(mesh_group)( + mesh_group.order, new_indices[group_nr], + new_nodes[group_nr], + unit_nodes=mesh_group.unit_nodes)) from meshmode.mesh import Mesh - part_mesh = Mesh(new_vertices, new_mesh_groups) + part_mesh = Mesh( + new_vertices, + new_mesh_groups, + is_conforming=mesh.is_conforming) return part_mesh, queried_elems @@ -318,7 +323,10 @@ def perform_flips(mesh, flip_flags, skip_tests=False): new_groups.append(new_grp) - return Mesh(mesh.vertices, new_groups, skip_tests=skip_tests) + return Mesh( + mesh.vertices, new_groups, skip_tests=skip_tests, + is_conforming=mesh.is_conforming, + ) # }}} @@ -435,7 +443,10 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False): from meshmode.mesh import Mesh return Mesh(vertices, new_groups, skip_tests=skip_tests, nodal_adjacency=nodal_adjacency, - facial_adjacency_groups=facial_adjacency_groups) + facial_adjacency_groups=facial_adjacency_groups, + is_conforming=all( + mesh.is_conforming + for mesh in meshes)) # }}} @@ -467,7 +478,8 @@ def map_mesh(mesh, f): # noqa from meshmode.mesh import Mesh return Mesh(vertices, new_groups, skip_tests=True, nodal_adjacency=mesh.nodal_adjacency_init_arg(), - facial_adjacency_groups=mesh._facial_adjacency_groups) + facial_adjacency_groups=mesh._facial_adjacency_groups, + is_conforming=mesh.is_conforming) # }}} diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index 9f3529f8..28638543 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -65,24 +65,34 @@ class TreeRayNode(object): self.adjacent_add_diff = [] -class Refiner(object): +class _Tesselation(RecordWithoutPickling): + + def __init__(self, children, ref_vertices): + RecordWithoutPickling.__init__(self, + children=children, + ref_vertices=ref_vertices,) - class _Tesselation(RecordWithoutPickling): - def __init__(self, children, ref_vertices): - RecordWithoutPickling.__init__(self, - ref_vertices=ref_vertices, children=children) +class _GroupRefinementRecord(RecordWithoutPickling): - class _GroupRefinementRecord(RecordWithoutPickling): + def __init__(self, tesselation, element_mapping): + RecordWithoutPickling.__init__(self, + tesselation=tesselation, element_mapping=element_mapping) - def __init__(self, tesselation, element_mapping): - RecordWithoutPickling.__init__(self, - tesselation=tesselation, element_mapping=element_mapping) + +class Refiner(object): # {{{ constructor def __init__(self, mesh): - from meshmode.mesh.tesselate import tesselateseg, tesselatetet, tesselatetri + if mesh.is_conforming is not True: + raise ValueError("Refiner can only be used with meshes that are known " + "to be conforming. If you would like to refine non-conforming " + "meshes and do not need adjacency information, consider " + "using RefinerWithoutAdjacency.") + + from meshmode.mesh.refinement.tesselate import \ + tesselateseg, tesselatetet, tesselatetri self.lazy = False self.seen_tuple = {} self.group_refinement_records = [] @@ -591,7 +601,7 @@ class Refiner(object): from meshmode.mesh.refinement.resampler import ( SimplexResampler) resampler = SimplexResampler() - tesselation = self._Tesselation( + tesselation = _Tesselation( self.simplex_result[grp.dim], self.simplex_node_tuples[grp.dim]) else: @@ -678,7 +688,7 @@ class Refiner(object): # if len(cur_list[len(cur_list)-1]) self.group_refinement_records.append( - self._GroupRefinementRecord(tesselation, element_mapping)) + _GroupRefinementRecord(tesselation, element_mapping)) #clear connectivity data for grp in self.last_mesh.groups: @@ -766,13 +776,18 @@ class Refiner(object): from meshmode.mesh import Mesh + refine_flags = refine_flags.astype(np.bool) + self.previous_mesh = self.last_mesh self.last_mesh = Mesh( vertices, new_mesh_el_groups, nodal_adjacency=self.generate_nodal_adjacency( totalnelements, nvertices, groups), vertex_id_dtype=self.last_mesh.vertex_id_dtype, - element_id_dtype=self.last_mesh.element_id_dtype) + element_id_dtype=self.last_mesh.element_id_dtype, + is_conforming=( + self.last_mesh.is_conforming + and (refine_flags.all() or (~refine_flags).all()))) return self.last_mesh # }}} diff --git a/meshmode/mesh/tesselate.py b/meshmode/mesh/refinement/tesselate.py similarity index 57% rename from meshmode/mesh/tesselate.py rename to meshmode/mesh/refinement/tesselate.py index 6b2315ed..ed80248c 100644 --- a/meshmode/mesh/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -1,3 +1,31 @@ +from __future__ import division, print_function + +__copyright__ = """ +Copyright (C) 2018 Andreas Kloeckner +Copyright (C) 2014-6 Shivam Gupta +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ as gnitstam @@ -6,10 +34,20 @@ def add_tuples(a, b): return tuple(ac+bc for ac, bc in zip(a, b)) +def halve_tuple(a): + def halve(x): + d, r = divmod(x, 2) + if r: + raise ValueError("%s is not evenly divisible by two" % x) + return d + + return tuple(halve(ac) for ac in a) + + def tesselateseg(): node_tuples = [(0,), (1,), (2,)] result = [(0, 1), (1, 2)] - return [node_tuples, result] + return node_tuples, result def tesselatetri(): @@ -39,7 +77,7 @@ def tesselatetri(): # positively oriented try_add_tri(current, (0, 0), (1, 0), (0, 1)) try_add_tri(current, (1, 0), (1, 1), (0, 1)) - return [node_tuples, result] + return node_tuples, result def tesselatetet(): @@ -77,4 +115,15 @@ def tesselatetet(): try_add_tet(current, (0, 1, 1), (0, 1, 0), (1, 1, 0), (1, 0, 1)) try_add_tet(current, (0, 1, 1), (1, 1, 1), (1, 0, 1), (1, 1, 0)) - return [node_tuples, result] + return node_tuples, result + + +def tesselate_simplex_bisection(dim): + if dim == 1: + return tesselateseg() + elif dim == 2: + return tesselatetri() + elif dim == 3: + return tesselatetet() + else: + raise ValueError("cannot tesselate %d-simplex" % dim) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 7009ae77..3ee4953f 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -590,7 +590,7 @@ def test_sanity_single_element(ctx_getter, dim, order, visualize=False): nodes=mp.warp_and_blend_nodes(dim, order).reshape(dim, 1, -1), dim=dim) - mesh = Mesh(vertices, [mg], nodal_adjacency=None, facial_adjacency_groups=None) + mesh = Mesh(vertices, [mg], is_conforming=True) from meshmode.discretization import Discretization from meshmode.discretization.poly_element import \ @@ -1063,10 +1063,11 @@ def test_ChainedDiscretizationConnection(ctx_getter): # noqa connections = [] + from meshmode.mesh.refinement import Refiner + refiner = Refiner(mesh) + def refine_discr(discr): mesh = discr.mesh - from meshmode.mesh.refinement import Refiner - refiner = Refiner(mesh) flags = refine_flags(mesh) refiner.refine(flags) connections.append( diff --git a/test/test_refinement.py b/test/test_refinement.py index ad47575c..d17d3db2 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -76,6 +76,10 @@ def even_refine_flags(spacing, mesh): return flags +def empty_refine_flags(mesh): + return np.zeros(mesh.nelements) + + def uniform_refine_flags(mesh): return np.ones(mesh.nelements) -- GitLab