diff --git a/.gitignore b/.gitignore index e090b1d87d5c4921b7bd12087cccd1b157c9cf07..1aed2289a06372e47ff275d875fa5819142d3529 100644 --- a/.gitignore +++ b/.gitignore @@ -15,3 +15,4 @@ a.out *.vtu .cache +.pytest_cache diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 19ac48c22df2c1b75acbf804fa89d966d4525690..2cf1358040411c548219fd85fdcef60f58a34d54 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -41,30 +41,30 @@ def _build_interpolation_batches_for_group( into is determined by where the refined unit nodes live relative to the coarse reference element. - For instance, consider the following refinement: + For instance, consider the following refinement:: - ______ ______ - |\ | |\ e| - | \ | |d\ | - | \ | |__\ | - | \ | => |\c|\ | - | \ | |a\|b\ | - | \| | | \| - ‾‾‾‾‾‾ ‾‾‾‾‾‾ + ______ ______ + |\ | |\ e| + | \ | |d\ | + | \ | |__\ | + | \ | => |\c|\ | + | \ | |a\|b\ | + | \| | | \| + ‾‾‾‾‾‾ ‾‾‾‾‾‾ Here, the discretization unit nodes for elements a,b,c,d,e will each have different positions relative to the reference element, so each element gets its own batch. On the other - hand, for + hand, for:: - ______ ______ - |\ | |\ f|\e| - | \ | |d\ |g\| - | \ | |__\|__| - | \ | => |\c|\ | - | \ | |a\|b\h| - | \| | | \| - ‾‾‾‾‾‾ ‾‾‾‾‾‾ + ______ ______ + |\ | |\ f|\e| + | \ | |d\ |g\| + | \ | |__\|__| + | \ | => |\c|\ | + | \ | |a\|b\h| + | \| | | \| + ‾‾‾‾‾‾ ‾‾‾‾‾‾ the pairs {a,e}, {b,f}, {c,g}, {d,h} can share interpolation batches because their unit nodes are mapped from the same part @@ -131,10 +131,11 @@ def make_refinement_connection(refiner, coarse_discr, group_factory): DirectDiscretizationConnection) coarse_mesh = refiner.get_previous_mesh() - fine_mesh = refiner.last_mesh + fine_mesh = refiner.get_current_mesh() + if coarse_discr.mesh != coarse_mesh: raise ValueError( - "course_discr must live on the same mesh given to the refiner!") + "coarse_discr does not live on the same mesh given to the refiner") from meshmode.discretization import Discretization fine_discr = Discretization( diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index d3368b440ea44a862f8dc694f986e69455e65d2a..89e97a137fd84b1ed291d17f80192afb6854b09e 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 198c45f4de53923153a28246be16cdd836cfd3f4..ac4928b748d1ed0b1d1cc265fdc49d31ee80240c 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,15 +560,16 @@ 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) - from meshmode.mesh.refinement import Refiner + from meshmode.mesh.refinement import RefinerWithoutAdjacency - refiner = Refiner(unwarped_mesh) + # These come out conformal, so we're OK to use the faster refiner. + refiner = RefinerWithoutAdjacency(unwarped_mesh) for i in range(uniform_refinement_rounds): refiner.refine_uniformly() @@ -730,8 +727,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) # }}} @@ -805,11 +801,11 @@ def warp_and_refine_until_resolved( from modepy.modes import simplex_onb from modepy.matrices import vandermonde from modepy.modal_decay import simplex_interp_error_coefficient_estimator_matrix - from meshmode.mesh.refinement import Refiner + from meshmode.mesh.refinement import Refiner, RefinerWithoutAdjacency logger.info("warp_and_refine_until_resolved: start") - if isinstance(unwarped_mesh_or_refiner, Refiner): + if isinstance(unwarped_mesh_or_refiner, (Refiner, RefinerWithoutAdjacency)): refiner = unwarped_mesh_or_refiner unwarped_mesh = refiner.get_current_mesh() else: diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index fd8a82a288ca4536cb9f65d9c3fad617c319a82c..c663a3eb6ac863aceb89b4f5d4ed1ea5f572a246 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 664355e0f043089af35efbd1ed79b0fdec7e71d8..36c38eda9d6704f5f59a6e589578d019d7485b88 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 9f3529f8cf4a79a1360fd6cb98642278eb4551ff..0d7b29b31d57387063a907e14f8f961a02e05d7a 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -27,6 +27,8 @@ import numpy as np import itertools from six.moves import range from pytools import RecordWithoutPickling +from meshmode.mesh.refinement.no_adjacency import ( # noqa: F401 + RefinerWithoutAdjacency) import logging @@ -34,6 +36,7 @@ logger = logging.getLogger(__name__) __doc__ = """ .. autoclass :: Refiner +.. autoclass :: RefinerWithoutAdjacency .. autofunction :: refine_uniformly """ @@ -65,24 +68,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): +class _GroupRefinementRecord(RecordWithoutPickling): - def __init__(self, children, ref_vertices): - RecordWithoutPickling.__init__(self, - ref_vertices=ref_vertices, children=children) + def __init__(self, tesselation, element_mapping): + RecordWithoutPickling.__init__(self, + tesselation=tesselation, element_mapping=element_mapping) - class _GroupRefinementRecord(RecordWithoutPickling): - 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 +604,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 +691,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 +779,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 # }}} @@ -926,21 +944,4 @@ def refine_uniformly(mesh, iterations): return mesh -class WarpingRefinerWrapper(object): - """A refiner that wraps an :attr:`inner_refiner` operating on an - un-warped mesh, that, for all external mesh retrieval purposes, - returns the warped mesh. - """ - - def __init__(self, inner_refiner, warp_mesh): - self.inner_refiner = inner_refiner - self.warp_mesh = warp_mesh - - def refine(self, refine_flags): - mesh = self.inner_refiner.refine(refine_flags) - return self.warp_mesh(mesh) - - def get_current_mesh(self): - return self.warp_mesh(self.inner_refiner.get_current_mesh()) - # vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py new file mode 100644 index 0000000000000000000000000000000000000000..4413117d0ea69aef016c353b24a16011cc55c9be --- /dev/null +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -0,0 +1,306 @@ +from __future__ import division, print_function + +__copyright__ = """ +Copyright (C) 2018 Andreas Kloeckner +Copyright (C) 2014-6 Shivam Gupta +Copyright (C) 2016 Matt Wala +""" + +__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. +""" + + +import numpy as np +from six.moves import range # noqa: F401 +from pytools import RecordWithoutPickling + +from pytools import memoize_method + + +import logging +logger = logging.getLogger(__name__) + + +class _TesselationInfo(RecordWithoutPickling): + + def __init__(self, children, ref_vertices, orig_vertex_indices, + midpoint_indices, midpoint_vertex_pairs, resampler): + RecordWithoutPickling.__init__(self, + children=children, + ref_vertices=ref_vertices, + orig_vertex_indices=orig_vertex_indices, + midpoint_indices=midpoint_indices, + midpoint_vertex_pairs=midpoint_vertex_pairs, + resampler=resampler) + + +class _GroupRefinementRecord(RecordWithoutPickling): + + def __init__(self, tesselation, element_mapping): + RecordWithoutPickling.__init__(self, + tesselation=tesselation, element_mapping=element_mapping) + + +class RefinerWithoutAdjacency(object): + """A refiner that may be applied to non-conforming + :class:`meshmode.mesh.Mesh` instances. It does not generate adjacency + information, and it is typically faster than + :class:`meshmode.mesh.refinement.Refiner`. + + .. note:: + + If the input meshes to this refiner are not conforming, then + the resulting meshes may contain duplicated vertices. + (I.e. two different numbers referring to the same geometric + vertex.) + """ + + def __init__(self, mesh): + self._current_mesh = mesh + self._previous_mesh = None + self.group_refinement_records = None + self.global_vertex_pair_to_midpoint = {} + + # {{{ build tesselation info + + @memoize_method + def _get_bisection_tesselation_info(self, group_type, dim): + from meshmode.mesh import SimplexElementGroup + if issubclass(group_type, SimplexElementGroup): + from meshmode.mesh.refinement.tesselate import \ + tesselate_simplex_bisection, add_tuples, halve_tuple + ref_vertices, children = tesselate_simplex_bisection(dim) + + orig_vertex_tuples = [(0,) * dim] + [ + (0,) * i + (2,) + (0,) * (dim-i-1) + for i in range(dim)] + node_dict = dict( + (ituple, idx) + for idx, ituple in enumerate(ref_vertices)) + orig_vertex_indices = [node_dict[vt] for vt in orig_vertex_tuples] + + from meshmode.mesh.refinement.resampler import SimplexResampler + resampler = SimplexResampler() + vertex_pair_to_midpoint_order = \ + resampler.get_vertex_pair_to_midpoint_order(dim) + + midpoint_idx_to_vertex_pair = {} + for vpair, mpoint_idx in vertex_pair_to_midpoint_order.items(): + midpoint_idx_to_vertex_pair[mpoint_idx] = vpair + + midpoint_vertex_pairs = [ + midpoint_idx_to_vertex_pair[i] + for i in range(len(midpoint_idx_to_vertex_pair))] + + midpoint_indices = [ + node_dict[ + halve_tuple( + add_tuples( + orig_vertex_tuples[v1], + orig_vertex_tuples[v2]))] + for v1, v2 in midpoint_vertex_pairs] + + return _TesselationInfo( + ref_vertices=ref_vertices, + children=np.array(children), + orig_vertex_indices=np.array(orig_vertex_indices), + midpoint_indices=np.array(midpoint_indices), + midpoint_vertex_pairs=midpoint_vertex_pairs, + resampler=resampler, + ) + + else: + raise NotImplementedError( + "bisection for elements groups of type %s" + % group_type.__name__) + + # }}} + + def refine_uniformly(self): + flags = np.ones(self._last_mesh.nelements, dtype=bool) + self.refine(flags) + + # {{{ refinement top-level + + def refine(self, refine_flags): + """ + :arg refine_flags: a :class:`numpy.ndarray` of dtype bool of length + ``mesh.nelements`` indicating which elements should be split. + """ + + mesh = self._current_mesh + + refine_flags = np.asarray(refine_flags, dtype=np.bool) + + if len(refine_flags) != mesh.nelements: + raise ValueError("length of refine_flags does not match " + "element count of last generated mesh") + + new_el_groups = [] + group_refinement_records = [] + additional_vertices = [] + inew_vertex = mesh.nvertices + + for igrp, group in enumerate(mesh.groups): + bisection_info = self._get_bisection_tesselation_info( + type(group), group.dim) + + # {{{ compute counts and index arrays + + grp_flags = refine_flags[ + group.element_nr_base: + group.element_nr_base+group.nelements] + + nchildren = len(bisection_info.children) + nchild_elements = np.ones(group.nelements, dtype=mesh.element_id_dtype) + nchild_elements[grp_flags] = nchildren + + child_el_indices = np.empty( + group.nelements+1, dtype=mesh.element_id_dtype) + child_el_indices[0] = 0 + child_el_indices[1:] = np.cumsum(nchild_elements) + + unrefined_el_new_indices = child_el_indices[:-1][~grp_flags] + refining_el_old_indices, = np.where(grp_flags) + + new_nelements = child_el_indices[-1] + + # }}} + + group_refinement_records.append( + _GroupRefinementRecord( + tesselation=bisection_info, + element_mapping=[ + list(range( + child_el_indices[iel], + child_el_indices[iel]+nchild_elements[iel])) + for iel in range(group.nelements)])) + + # {{{ get new vertices together + + midpoints = bisection_info.resampler.get_midpoints( + group, bisection_info, refining_el_old_indices) + + new_vertex_indices = np.empty( + (new_nelements, group.vertex_indices.shape[1]), + dtype=mesh.vertex_id_dtype) + new_vertex_indices.fill(-17) + + # copy over unchanged vertices + new_vertex_indices[unrefined_el_new_indices] = \ + group.vertex_indices[~grp_flags] + + for old_iel in refining_el_old_indices: + new_iel_base = child_el_indices[old_iel] + + refining_vertices = np.empty(len(bisection_info.ref_vertices), + dtype=mesh.vertex_id_dtype) + refining_vertices.fill(-17) + + # carry over old vertices + refining_vertices[bisection_info.orig_vertex_indices] = \ + group.vertex_indices[old_iel] + + for imidpoint, (iref_midpoint, (v1, v2)) in enumerate(zip( + bisection_info.midpoint_indices, + bisection_info.midpoint_vertex_pairs)): + + global_v1 = group.vertex_indices[old_iel, v1] + global_v2 = group.vertex_indices[old_iel, v2] + + if global_v1 > global_v2: + global_v1, global_v2 = global_v2, global_v1 + + try: + global_midpoint = self.global_vertex_pair_to_midpoint[ + global_v1, global_v2] + except KeyError: + global_midpoint = inew_vertex + additional_vertices.append(midpoints[old_iel][:, imidpoint]) + inew_vertex += 1 + + refining_vertices[iref_midpoint] = global_midpoint + + assert (refining_vertices >= 0).all() + + new_vertex_indices[new_iel_base:new_iel_base+nchildren] = \ + refining_vertices[bisection_info.children] + + assert (new_vertex_indices >= 0).all() + + # }}} + + # {{{ get new nodes together + + new_nodes = np.empty( + (mesh.ambient_dim, new_nelements, group.nunit_nodes), + dtype=group.nodes.dtype) + + new_nodes.fill(float("nan")) + + # copy over unchanged nodes + new_nodes[:, unrefined_el_new_indices] = group.nodes[:, ~grp_flags] + + tesselated_nodes = bisection_info.resampler.get_tesselated_nodes( + group, bisection_info, refining_el_old_indices) + + for old_iel in refining_el_old_indices: + new_iel_base = child_el_indices[old_iel] + new_nodes[:, new_iel_base:new_iel_base+nchildren, :] = \ + tesselated_nodes[old_iel] + + assert (~np.isnan(new_nodes)).all() + + # }}} + + new_el_groups.append( + type(group)( + order=group.order, + vertex_indices=new_vertex_indices, + nodes=new_nodes, + unit_nodes=group.unit_nodes)) + + new_vertices = np.empty( + (mesh.ambient_dim, mesh.nvertices + len(additional_vertices)), + mesh.vertices.dtype) + new_vertices[:, :mesh.nvertices] = mesh.vertices + new_vertices[:, mesh.nvertices:] = np.array(additional_vertices).T + + from meshmode.mesh import Mesh + new_mesh = Mesh(new_vertices, new_el_groups, is_conforming=( + mesh.is_conforming + and (refine_flags.all() or (~refine_flags).all()))) + + self.group_refinement_records = group_refinement_records + self._current_mesh = new_mesh + self._previous_mesh = mesh + + return new_mesh + + # }}} + + def get_current_mesh(self): + return self._current_mesh + + def get_previous_mesh(self): + return self._previous_mesh + + +# vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/resampler.py b/meshmode/mesh/refinement/resampler.py index cea22b3aa3f7c36b23f823e06646bb27e2723110..794c1c5a2362d1255f5e3d329c48fb766c6229a9 100644 --- a/meshmode/mesh/refinement/resampler.py +++ b/meshmode/mesh/refinement/resampler.py @@ -43,7 +43,8 @@ class SimplexResampler(object): :func:`meshmode.mesh.tesselate.tesselatetet()`. """ - def get_vertex_pair_to_midpoint_order(self, dim): + @staticmethod + def get_vertex_pair_to_midpoint_order(dim): """ :arg dim: Dimension of the element @@ -58,7 +59,8 @@ class SimplexResampler(object): range(nmidpoints) )) - def get_midpoints(self, group, tesselation, elements): + @staticmethod + def get_midpoints(group, tesselation, elements): """ Compute the midpoints of the vertices of the specified elements. @@ -89,7 +91,8 @@ class SimplexResampler(object): return dict(zip(elements, resamp_midpoints)) - def get_tesselated_nodes(self, group, tesselation, elements): + @staticmethod + def get_tesselated_nodes(group, tesselation, elements): """ Compute the nodes of the child elements according to the tesselation. 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 6b2315ed053eb3ff37a02dd74fe8c02a5eb4aefc..ed80248cdbfffa34004c8bc218f0aa47f984f02d 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 7009ae77348c56388f15218535ba5f4937c0fcb6..3ee4953f2d2876c625ebc4819d08d22eb4a82ff2 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 ad47575ceeb32115baf0a44a0285b5308399bd38..8c8ceee02a4de6aabe92bc8c524272f79714dc56 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -33,7 +33,7 @@ from pyopencl.tools import ( # noqa from meshmode.mesh.generation import ( # noqa generate_icosahedron, generate_box_mesh, make_curve_mesh, ellipse) from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry -from meshmode.mesh.refinement import Refiner +from meshmode.mesh.refinement import Refiner, RefinerWithoutAdjacency from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, @@ -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) @@ -148,6 +152,10 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): check_nodal_adj_against_geometry(mesh) +@pytest.mark.parametrize("refiner_cls", [ + Refiner, + RefinerWithoutAdjacency + ]) @pytest.mark.parametrize("group_factory", [ InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, @@ -166,9 +174,10 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): #partial(random_refine_flags, 0.4) partial(even_refine_flags, 2) ]) +# test_refinement_connection(cl._csc, RefinerWithoutAdjacency, PolynomialWarpAndBlendGroupFactory, 'warp', 2, [4, 5, 6], 5, partial(even_refine_flags, 2)) # noqa: E501 def test_refinement_connection( - ctx_getter, group_factory, mesh_name, dim, mesh_pars, mesh_order, - refine_flags, plot_mesh=False): + ctx_getter, refiner_cls, group_factory, + mesh_name, dim, mesh_pars, mesh_order, refine_flags, visualize=False): from random import seed seed(13) @@ -215,7 +224,7 @@ def test_refinement_connection( discr = Discretization(cl_ctx, mesh, group_factory(order)) - refiner = Refiner(mesh) + refiner = refiner_cls(mesh) flags = refine_flags(mesh) refiner.refine(flags) @@ -231,7 +240,7 @@ def test_refinement_connection( f_interp = connection(queue, f_coarse).with_queue(queue) f_true = f(x_fine).with_queue(queue) - if plot_mesh: + if visualize == "dots": import matplotlib.pyplot as plt x = x.get(queue) err = np.array(np.log10( @@ -243,6 +252,16 @@ def test_refinement_connection( plt.colorbar(cmap) plt.show() + elif visualize == "vtk": + from meshmode.discretization.visualization import make_visualizer + fine_vis = make_visualizer(queue, fine_discr, mesh_order) + + fine_vis.write_vtk_file( + "refine-fine-%s-%dd-%s.vtu" % (mesh_name, dim, mesh_par), [ + ("f_interp", f_interp), + ("f_true", f_true), + ]) + import numpy.linalg as la err = la.norm((f_interp - f_true).get(queue), np.inf) eoc_rec.add_data_point(h, err)