From 1b1b434b507d23f684d799c6213bb777a16ee338 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 3 Dec 2020 23:48:02 -0600 Subject: [PATCH 01/20] start implementing refinement --- meshmode/mesh/refinement/__init__.py | 13 +- meshmode/mesh/refinement/no_adjacency.py | 218 ++++++++++++++-------- meshmode/mesh/refinement/resampler.py | 222 ++++++++++++++--------- meshmode/mesh/refinement/tesselate.py | 4 + meshmode/mesh/refinement/utils.py | 1 - test/test_refinement.py | 12 +- 6 files changed, 298 insertions(+), 172 deletions(-) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index be101b487..f501985a4 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -27,6 +27,7 @@ from meshmode.mesh.refinement.no_adjacency import ( # noqa: F401 RefinerWithoutAdjacency) +from meshmode.mesh.refinement.resampler import TesselationInfo as _Tesselation import logging logger = logging.getLogger(__name__) @@ -65,14 +66,6 @@ def __init__(self, left_vertex, right_vertex, adjacent_elements=[]): self.adjacent_add_diff = [] -class _Tesselation(RecordWithoutPickling): - - def __init__(self, children, ref_vertices): - RecordWithoutPickling.__init__(self, - children=children, - ref_vertices=ref_vertices,) - - class _GroupRefinementRecord(RecordWithoutPickling): def __init__(self, tesselation, element_mapping): @@ -617,8 +610,8 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr SimplexResampler) resampler = SimplexResampler() tesselation = _Tesselation( - self.simplex_result[grp.dim], - self.simplex_node_tuples[grp.dim]) + children=self.simplex_result[grp.dim], + ref_vertices=self.simplex_node_tuples[grp.dim]) else: raise NotImplementedError("unimplemented: midpoint finding" "for non simplex elements") diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index adf1fc2f9..9a94ac0a2 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -26,33 +26,156 @@ import numpy as np -from pytools import RecordWithoutPickling - -from pytools import memoize_method +import modepy as mp +from pytools import RecordWithoutPickling +from pytools import memoize, memoize_method import logging logger = logging.getLogger(__name__) +# {{{ tesselation + class _TesselationInfo(RecordWithoutPickling): + """ + .. attribute:: shape + + .. attribute:: children + .. attribute:: ref_vertices - 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) + .. attribute:: orig_vertex_indices + .. attribute:: midpoint_indices + .. attribute:: midpoint_vertex_pairs + """ class _GroupRefinementRecord(RecordWithoutPickling): + """ + .. attribute:: tesselation + .. attribute:: element_mapping + """ + - def __init__(self, tesselation, element_mapping): - RecordWithoutPickling.__init__(self, - tesselation=tesselation, element_mapping=element_mapping) +def _old_simplex_tesselation_info(dim): + 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 = { + 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, + ) + + +def _midpoint_tuples(a, b): + def halve(ai, bi): + d, r = divmod(ai + bi, 2) + if r: + raise ValueError(f"({ai} + {bi}) is not evenly divisible by two") + + return d + + return tuple(halve(ai, bi) for ai, bi in zip(a, b)) + + +def get_simplex_tesselation_info(shape): + # ref_vertices refer to the nodes in the 2nd order simplex + # + # F + # | \ + # | \ + # D----E + # | /| \ + # | / | \ + # A----B----C + # + # orig_vertices is just (A, C, F) + # + # we then find the midpoints (B, D, E) and the vertex pairs around them. + + space = mp.space_for_shape(shape, 2) + ref_vertices = mp.node_tuples_for_space(space) + ref_vertices_to_index = {rv: i for i, rv in enumerate(ref_vertices)} + + from pytools import add_tuples + space = mp.space_for_shape(shape, 1) + orig_vertices = tuple([ + add_tuples(vt, vt) for vt in mp.node_tuples_for_space(space) + ]) + orig_vertex_indices = [ref_vertices_to_index[vt] for vt in orig_vertices] + + from meshmode.mesh.refinement.resampler import get_ref_midpoints + midpoints = get_ref_midpoints(shape, ref_vertices) + midpoint_indices = [ref_vertices_to_index[mp] for mp in midpoints] + + midpoint_to_vertex_pairs = { + midpoint: (i, j) + for i, ivt in enumerate(orig_vertices) + for j, jvt in enumerate(orig_vertices) + for midpoint in [_midpoint_tuples(ivt, jvt)] + if i < j and midpoint in midpoints + } + midpoint_vertex_pairs = [midpoint_to_vertex_pairs[m] for m in midpoints] + + return _TesselationInfo( + shape=shape, + ref_vertices=ref_vertices, + children=np.array(mp.submesh_for_shape(shape, ref_vertices)), + orig_vertex_indices=np.array(orig_vertex_indices), + midpoint_indices=np.array(midpoint_indices), + midpoint_vertex_pairs=midpoint_vertex_pairs, + ) + + +def get_hypercube_tesselation_info(shape): + return get_simplex_tesselation_info(shape) + +@memoize +def get_bisection_tesselation_info(group_type, dim): + from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup + if issubclass(group_type, SimplexElementGroup): + return get_simplex_tesselation_info(mp.Simplex(dim)) + elif issubclass(group_type, TensorProductElementGroup): + return get_hypercube_tesselation_info(mp.Hypercube(dim)) + else: + raise NotImplementedError( + "bisection for element groups of type {group_type.__name__}") + +# }}} class RefinerWithoutAdjacency: @@ -81,61 +204,6 @@ def __init__(self, mesh): 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 = { - 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._current_mesh.nelements, dtype=bool) return self.refine(flags) @@ -163,9 +231,11 @@ def refine(self, refine_flags): if perform_vertex_updates: inew_vertex = mesh.nvertices + from meshmode.mesh.refinement.resampler import \ + get_midpoints, get_tesselated_nodes + for igrp, group in enumerate(mesh.groups): - bisection_info = self._get_bisection_tesselation_info( - type(group), group.dim) + bisection_info = get_bisection_tesselation_info(type(group), group.dim) # {{{ compute counts and index arrays @@ -201,7 +271,7 @@ def refine(self, refine_flags): # {{{ get new vertices together if perform_vertex_updates: - midpoints = bisection_info.resampler.get_midpoints( + midpoints = get_midpoints(bisection_info.shape, group, bisection_info, refining_el_old_indices) new_vertex_indices = np.empty( @@ -269,7 +339,7 @@ def refine(self, refine_flags): # copy over unchanged nodes new_nodes[:, unrefined_el_new_indices] = group.nodes[:, ~grp_flags] - tesselated_nodes = bisection_info.resampler.get_tesselated_nodes( + tesselated_nodes = get_tesselated_nodes(bisection_info.shape, group, bisection_info, refining_el_old_indices) for old_iel in refining_el_old_indices: diff --git a/meshmode/mesh/refinement/resampler.py b/meshmode/mesh/refinement/resampler.py index fd6ae9396..39cd539df 100644 --- a/meshmode/mesh/refinement/resampler.py +++ b/meshmode/mesh/refinement/resampler.py @@ -20,9 +20,11 @@ THE SOFTWARE. """ +from functools import singledispatch import numpy as np import modepy as mp +from pytools import RecordWithoutPickling import logging logger = logging.getLogger(__name__) @@ -30,108 +32,160 @@ # {{{ resampling simplex points for refinement -# NOTE: Class internal to refiner: do not make documentation public. class SimplexResampler: - """ - Resampling of points on simplex elements for refinement. - - Most methods take a ``tesselation`` parameter. - The tesselation should follow the format of - :func:`meshmode.mesh.tesselate.tesselatetri()` or - :func:`meshmode.mesh.tesselate.tesselatetet()`. - """ - @staticmethod def get_vertex_pair_to_midpoint_order(dim): - """ - :arg dim: Dimension of the element - - :return: A :class:`dict` mapping the vertex pair :math:`(v1, v2)` (with - :math:`v1 < v2`) to the number of the midpoint in the tesselation - ordering (the numbering is restricted to the midpoints, so there - are no gaps in the numbering) - """ - nmidpoints = dim * (dim + 1) // 2 - return dict(zip( - ((i, j) for j in range(dim + 1) for i in range(j)), - range(nmidpoints) - )) + return get_vertex_pair_to_midpoint_order(mp.Simplex(dim)) @staticmethod def get_midpoints(group, tesselation, elements): - """ - Compute the midpoints of the vertices of the specified elements. - - :arg group: An instance of :class:`meshmode.mesh.SimplexElementGroup` - :arg tesselation: With attributes `ref_vertices`, `children` - :arg elements: A list of (group-relative) element numbers - - :return: A :class:`dict` mapping element numbers to midpoint - coordinates, with each value in the map having shape - ``(ambient_dim, nmidpoints)``. The ordering of the midpoints - follows their ordering in the tesselation (see also - :meth:`SimplexResampler.get_vertex_pair_to_midpoint_order`) - """ - if group.vertex_indices is not None: - assert len(group.vertex_indices[0]) == group.dim + 1 - - # Get midpoints, converted to unit coordinates. - midpoints = -1 + np.array([vertex for vertex in - tesselation.ref_vertices if 1 in vertex], dtype=float) - - resamp_mat = mp.resampling_matrix( - mp.simplex_best_available_basis(group.dim, group.order), + return get_midpoints(mp.Simplex(group.dim), + group, tesselation, elements) + + @staticmethod + def get_tesselated_nodes(group, tesselation, elements): + return get_tesselated_nodes(mp.Simplex(group.dim), + group, tesselation, elements) + +# }}} + + +# {{{ interface + +# NOTE: internal to refiners: do not make documentation public. + +class TesselationInfo(RecordWithoutPickling): + """ + .. attribute:: ref_vertices + .. attribute:: children + """ + +@singledispatch +def map_unit_nodes_to_children(shape: mp.Shape, tesselation, unit_nodes): + raise NotImplementedError + + +@singledispatch +def get_ref_midpoints(shape: mp.Shape, ref_vertices): + from pytools import add_tuples + space = mp.space_for_shape(shape, 1) + orig_vertices = [ + add_tuples(vt, vt) for vt in mp.node_tuples_for_space(space) + ] + return [rv for rv in ref_vertices if rv not in orig_vertices] + + +@singledispatch +def get_midpoints(shape: mp.Shape, group, tesselation, elements): + """Compute the midpoints of the vertices of the specified elements. + + :arg group: an instance of :class:`meshmode.mesh.MeshElementGroup`. + :arg tesselation: a :class:`TesselationInfo`. + :arg elements: a list of (group-relative) element numbers. + + :return: A :class:`dict` mapping element numbers to midpoint + coordinates, with each value in the map having shape + ``(ambient_dim, nmidpoints)``. The ordering of the midpoints + follows their ordering in the tesselation (see also + :meth:`SimplexResampler.get_vertex_pair_to_midpoint_order`) + """ + if shape.dim != group.dim: + raise ValueError("shape and group dimension do not match") + + if group.vertex_indices is not None: + assert len(group.vertex_indices[0]) == shape.nvertices + + # get midpoints, converted to unit coordinates. + midpoints = -1 + np.array(get_ref_midpoints(shape, tesselation.ref_vertices)) + + space = mp.space_for_shape(shape, group.order) + resampling_mat = mp.resampling_matrix( + mp.basis_for_space(space, shape).functions, midpoints.T, group.unit_nodes) - resamp_midpoints = np.einsum("mu,deu->edm", - resamp_mat, - group.nodes[:, elements]) + resampled_midpoints = np.einsum("mu,deu->edm", + resampling_mat, group.nodes[:, elements]) - return dict(zip(elements, resamp_midpoints)) + return dict(zip(elements, resampled_midpoints)) - @staticmethod - def get_tesselated_nodes(group, tesselation, elements): - """ - Compute the nodes of the child elements according to the tesselation. - - :arg group: An instance of :class:`meshmode.mesh.SimplexElementGroup` - :arg tesselation: With attributes `ref_vertices`, `children` - :arg elements: A list of (group-relative) element numbers - - :return: A :class:`dict` mapping element numbers to node - coordinates, with each value in the map having shape - ``(ambient_dim, nchildren, nunit_nodes)``. - The ordering of the child nodes follows the ordering - of ``tesselation.children.`` - """ - if group.vertex_indices is not None: - assert len(group.vertex_indices[0]) == group.dim + 1 - - from meshmode.mesh.refinement.utils import map_unit_nodes_to_children - - # Get child unit node coordinates. - child_unit_nodes = np.hstack(list( - map_unit_nodes_to_children(group.unit_nodes, tesselation))) - - resamp_mat = mp.resampling_matrix( - mp.simplex_best_available_basis(group.dim, group.order), + +@singledispatch +def get_tesselated_nodes(shape: mp.Shape, group, tesselation, elements): + """Compute the nodes of the child elements according to the tesselation. + + :arg group: An instance of :class:`meshmode.mesh.MeshElementGroup`. + :arg tesselation: a :class:`TesselationInfo`. + :arg elements: A list of (group-relative) element numbers. + + :return: A :class:`dict` mapping element numbers to node + coordinates, with each value in the map having shape + ``(ambient_dim, nchildren, nunit_nodes)``. + The ordering of the child nodes follows the ordering + of ``tesselation.children.`` + """ + if shape.dim != group.dim: + raise ValueError("shape and group dimension do not match") + + if group.vertex_indices is not None: + assert len(group.vertex_indices[0]) == shape.nvertices + + # get child unit node coordinates. + child_unit_nodes = np.hstack(list( + map_unit_nodes_to_children(shape, tesselation, group.unit_nodes) + )) + + space = mp.space_for_shape(shape, group.order) + resampling_mat = mp.resampling_matrix( + mp.basis_for_space(space, shape).functions, child_unit_nodes, group.unit_nodes) - resamp_unit_nodes = np.einsum("cu,deu->edc", - resamp_mat, - group.nodes[:, elements]) + resampled_unit_nodes = np.einsum("cu,deu->edc", + resampling_mat, group.nodes[:, elements]) + + ambient_dim = len(group.nodes) + nunit_nodes = len(group.unit_nodes[0]) - ambient_dim = len(group.nodes) - nunit_nodes = len(group.unit_nodes[0]) + return { + el: resampled_unit_nodes[iel].reshape((ambient_dim, -1, nunit_nodes)) + for iel, el in enumerate(elements) + } + raise NotImplementedError(type(shape).__name__) - return {elem: - resamp_unit_nodes[ielem].reshape( - (ambient_dim, -1, nunit_nodes)) - for ielem, elem in enumerate(elements)} +# }}} + + +# {{{ simplex + +@map_unit_nodes_to_children.register(mp.Simplex) +def _(shape: mp.Simplex, tesselation, unit_nodes): + ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float) + assert len(unit_nodes.shape) == 2 + + for child_element in tesselation.children: + origin = np.vstack(ref_vertices[child_element[0]]) + basis = ref_vertices.T[:, child_element[1:]] - origin + + yield basis.dot((unit_nodes + 1) / 2) + origin - 1 # }}} +# {{{ hypercube + +@map_unit_nodes_to_children.register(mp.Hypercube) +def _(shape: mp.Hypercube, tesselation, unit_nodes): + ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float) + assert len(unit_nodes.shape) == 2 + + basis_indices = [1, 2, 4][:shape.dim] + for child_element in tesselation.children: + origin = np.vstack(ref_vertices[child_element[0]]) + basis = ref_vertices.T[:, np.array(child_element)[indices]] - origin + + yield basis.dot((unit_nodes + 1) / 2) + origin - 1 + +# }}} + # vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index 211a87f08..b099fd261 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -28,6 +28,10 @@ as gnitstam +def mul_tuples(t, m): + return tuple(m * a for a in t) + + def add_tuples(a, b): return tuple(ac+bc for ac, bc in zip(a, b)) diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index be64eb07c..1b717a803 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -47,7 +47,6 @@ def map_unit_nodes_to_children(unit_nodes, tesselation): :arg tesselation: With attributes `ref_vertices`, `children` """ ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float) - assert len(unit_nodes.shape) == 2 for child_element in tesselation.children: diff --git a/test/test_refinement.py b/test/test_refinement.py index afbda13ca..9314c93fa 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -36,6 +36,7 @@ from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry from meshmode.mesh.refinement import Refiner, RefinerWithoutAdjacency +from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup from meshmode.discretization.poly_element import ( InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, @@ -289,12 +290,17 @@ def f(x): or eoc_rec.max_error() < 1e-14) -@pytest.mark.parametrize("with_adjacency", [True, False]) -def test_uniform_refinement(with_adjacency): +@pytest.mark.parametrize(("group_cls", "with_adjacency"), [ + (SimplexElementGroup, True), + (SimplexElementGroup, False), + (TensorProductElementGroup, False) + ]) +def test_uniform_refinement(group_cls, with_adjacency): make_mesh = partial(generate_box_mesh, ( np.linspace(0.0, 1.0, 2), np.linspace(0.0, 1.0, 3), - np.linspace(0.0, 1.0, 2)), order=4) + np.linspace(0.0, 1.0, 2)), + order=4, group_cls=group_cls) mesh = make_mesh() from meshmode.mesh.refinement import refine_uniformly From 028b37b8cd8b7e4e8edae9eaf2eca40c9eea5147 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 6 Dec 2020 11:07:04 -0600 Subject: [PATCH 02/20] move all refinement helpers to one file --- .../discretization/connection/refinement.py | 13 +- meshmode/mesh/__init__.py | 6 +- meshmode/mesh/refinement/__init__.py | 69 +++- meshmode/mesh/refinement/no_adjacency.py | 180 +--------- meshmode/mesh/refinement/resampler.py | 191 ----------- meshmode/mesh/refinement/tesselate.py | 319 ++++++++++++------ meshmode/mesh/refinement/utils.py | 160 --------- 7 files changed, 309 insertions(+), 629 deletions(-) delete mode 100644 meshmode/mesh/refinement/resampler.py delete mode 100644 meshmode/mesh/refinement/utils.py diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 1124e5ad9..f09683050 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -87,9 +87,16 @@ def _build_interpolation_batches_for_group( to_bin.append(child_idx) fine_unit_nodes = fine_discr_group.unit_nodes - from meshmode.mesh.refinement.utils import map_unit_nodes_to_children - mapped_unit_nodes = map_unit_nodes_to_children( - fine_unit_nodes, record.tesselation) + meg = fine_discr_group.mesh_el_group + + from meshmode.mesh import _ModepyElementGroup + if isinstance(meg, _ModepyElementGroup): + from meshmode.mesh.refinement.tesselate import map_unit_nodes_to_children + mapped_unit_nodes = map_unit_nodes_to_children( + meg._modepy_shape, fine_unit_nodes, record.tesselation) + else: + raise TypeError("cannot construct refinement connection for groups " + f"of type '{type(meg).__name__}'") from itertools import chain for from_bin, to_bin, unit_nodes in zip( diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index ea298b185..bd400aa0c 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -914,10 +914,8 @@ def __eq__(self, other): and self.groups == other.groups and self.vertex_id_dtype == other.vertex_id_dtype and self.element_id_dtype == other.element_id_dtype - and (self._nodal_adjacency - == other._nodal_adjacency) - and (self._facial_adjacency_groups - == other._facial_adjacency_groups) + and self._nodal_adjacency == other._nodal_adjacency + and self._facial_adjacency_groups == other._facial_adjacency_groups and self.boundary_tags == other.boundary_tags and self.is_conforming == other.is_conforming) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index f501985a4..53ecbd162 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -27,7 +27,8 @@ from meshmode.mesh.refinement.no_adjacency import ( # noqa: F401 RefinerWithoutAdjacency) -from meshmode.mesh.refinement.resampler import TesselationInfo as _Tesselation +from meshmode.mesh.refinement.tesselate import \ + TesselationInfo, GroupRefinementRecord import logging logger = logging.getLogger(__name__) @@ -38,6 +39,49 @@ .. autofunction :: refine_uniformly """ +# {{{ deprecated + +class SimplexResampler: + @staticmethod + def get_vertex_pair_to_midpoint_order(dim): + nmidpoints = dim * (dim + 1) // 2 + return dict(zip( + ((i, j) for j in range(dim + 1) for i in range(j)), + range(nmidpoints) + )) + + @staticmethod + def get_midpoints(group, tesselation, elements): + import meshmode.mesh.refinement.tesselate as tess + return tess.get_midpoints(group, tesselation, elements) + + @staticmethod + def get_tesselated_nodes(group, tesselation, elements): + import meshmode.mesh.refinement.tesselate as tess + return tess.get_tesselated_nodes(group, tesselation, elements) + + +def tesselate_simplex(dim): + import modepy as mp + shape = mp.Simplex(dim) + space = mp.space_for_shape(shape, 2) + + node_tuples = mp.node_tuples_for_space(space) + return node_tuples, mp.submesh_for_shape(shape, node_tuples) + + +def tesselateseg(): + return tesselate_simplex(1) + + +def tesselatetri(): + return tesselate_simplex(2) + + +def tesselatetet(): + return tesselate_simplex(3) + +# }}} class TreeRayNode: """Describes a ray as a tree, this class represents each node in this tree @@ -66,13 +110,6 @@ def __init__(self, left_vertex, right_vertex, adjacent_elements=[]): self.adjacent_add_diff = [] -class _GroupRefinementRecord(RecordWithoutPickling): - - def __init__(self, tesselation, element_mapping): - RecordWithoutPickling.__init__(self, - tesselation=tesselation, element_mapping=element_mapping) - - class Refiner: """An older that mostly succeeds at preserving adjacency across non-conformal refinement. @@ -99,8 +136,6 @@ def __init__(self, mesh): "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,6 +626,10 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr del self.group_refinement_records[:] for grp_idx, grp in enumerate(self.last_mesh.groups): + if not isinstance(grp, SimplexElementGroup): + raise TypeError("refinement not supported for groups of type " + f"'{type(grp).__name__}'") + iel_base = grp.element_nr_base # List of lists mapping element number to new element number(s). element_mapping = [] @@ -606,10 +645,8 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr if len(grp.vertex_indices[iel_grp]) == grp.dim + 1: midpoints_to_find.append(iel_grp) if not resampler: - from meshmode.mesh.refinement.resampler import ( - SimplexResampler) resampler = SimplexResampler() - tesselation = _Tesselation( + tesselation = TesselationInfo( children=self.simplex_result[grp.dim], ref_vertices=self.simplex_node_tuples[grp.dim]) else: @@ -696,7 +733,10 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr # if len(cur_list[len(cur_list)-1]) self.group_refinement_records.append( - _GroupRefinementRecord(tesselation, element_mapping)) + GroupRefinementRecord( + tesselation=tesselation, + element_mapping=element_mapping) + ) #clear connectivity data for grp in self.last_mesh.groups: @@ -753,7 +793,6 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr if to_resample: # if simplex if is_simplex: - from meshmode.mesh.refinement.resampler import SimplexResampler resampler = SimplexResampler() new_nodes = resampler.get_tesselated_nodes( prev_group, refinement_record.tesselation, to_resample) diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index 9a94ac0a2..484363c72 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -24,160 +24,14 @@ THE SOFTWARE. """ - import numpy as np -import modepy as mp -from pytools import RecordWithoutPickling -from pytools import memoize, memoize_method +from meshmode.mesh.refinement.tesselate import GroupRefinementRecord import logging logger = logging.getLogger(__name__) -# {{{ tesselation - -class _TesselationInfo(RecordWithoutPickling): - """ - .. attribute:: shape - - .. attribute:: children - .. attribute:: ref_vertices - - .. attribute:: orig_vertex_indices - .. attribute:: midpoint_indices - .. attribute:: midpoint_vertex_pairs - """ - - -class _GroupRefinementRecord(RecordWithoutPickling): - """ - .. attribute:: tesselation - .. attribute:: element_mapping - """ - - -def _old_simplex_tesselation_info(dim): - 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 = { - 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, - ) - - -def _midpoint_tuples(a, b): - def halve(ai, bi): - d, r = divmod(ai + bi, 2) - if r: - raise ValueError(f"({ai} + {bi}) is not evenly divisible by two") - - return d - - return tuple(halve(ai, bi) for ai, bi in zip(a, b)) - - -def get_simplex_tesselation_info(shape): - # ref_vertices refer to the nodes in the 2nd order simplex - # - # F - # | \ - # | \ - # D----E - # | /| \ - # | / | \ - # A----B----C - # - # orig_vertices is just (A, C, F) - # - # we then find the midpoints (B, D, E) and the vertex pairs around them. - - space = mp.space_for_shape(shape, 2) - ref_vertices = mp.node_tuples_for_space(space) - ref_vertices_to_index = {rv: i for i, rv in enumerate(ref_vertices)} - - from pytools import add_tuples - space = mp.space_for_shape(shape, 1) - orig_vertices = tuple([ - add_tuples(vt, vt) for vt in mp.node_tuples_for_space(space) - ]) - orig_vertex_indices = [ref_vertices_to_index[vt] for vt in orig_vertices] - - from meshmode.mesh.refinement.resampler import get_ref_midpoints - midpoints = get_ref_midpoints(shape, ref_vertices) - midpoint_indices = [ref_vertices_to_index[mp] for mp in midpoints] - - midpoint_to_vertex_pairs = { - midpoint: (i, j) - for i, ivt in enumerate(orig_vertices) - for j, jvt in enumerate(orig_vertices) - for midpoint in [_midpoint_tuples(ivt, jvt)] - if i < j and midpoint in midpoints - } - midpoint_vertex_pairs = [midpoint_to_vertex_pairs[m] for m in midpoints] - - return _TesselationInfo( - shape=shape, - ref_vertices=ref_vertices, - children=np.array(mp.submesh_for_shape(shape, ref_vertices)), - orig_vertex_indices=np.array(orig_vertex_indices), - midpoint_indices=np.array(midpoint_indices), - midpoint_vertex_pairs=midpoint_vertex_pairs, - ) - - -def get_hypercube_tesselation_info(shape): - return get_simplex_tesselation_info(shape) - -@memoize -def get_bisection_tesselation_info(group_type, dim): - from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup - if issubclass(group_type, SimplexElementGroup): - return get_simplex_tesselation_info(mp.Simplex(dim)) - elif issubclass(group_type, TensorProductElementGroup): - return get_hypercube_tesselation_info(mp.Hypercube(dim)) - else: - raise NotImplementedError( - "bisection for element groups of type {group_type.__name__}") - -# }}} - - class RefinerWithoutAdjacency: """A refiner that may be applied to non-conforming :class:`meshmode.mesh.Mesh` instances. It does not generate adjacency @@ -231,11 +85,13 @@ def refine(self, refine_flags): if perform_vertex_updates: inew_vertex = mesh.nvertices - from meshmode.mesh.refinement.resampler import \ - get_midpoints, get_tesselated_nodes + from meshmode.mesh.refinement.tesselate import ( + get_bisection_tesselation_info, + get_midpoints, + get_tesselated_nodes) for igrp, group in enumerate(mesh.groups): - bisection_info = get_bisection_tesselation_info(type(group), group.dim) + tesselation = get_bisection_tesselation_info(type(group), group.dim) # {{{ compute counts and index arrays @@ -243,7 +99,7 @@ def refine(self, refine_flags): group.element_nr_base: group.element_nr_base+group.nelements] - nchildren = len(bisection_info.children) + nchildren = len(tesselation.children) nchild_elements = np.ones(group.nelements, dtype=mesh.element_id_dtype) nchild_elements[grp_flags] = nchildren @@ -260,8 +116,8 @@ def refine(self, refine_flags): # }}} group_refinement_records.append( - _GroupRefinementRecord( - tesselation=bisection_info, + GroupRefinementRecord( + tesselation=tesselation, element_mapping=[ list(range( child_el_indices[iel], @@ -271,8 +127,8 @@ def refine(self, refine_flags): # {{{ get new vertices together if perform_vertex_updates: - midpoints = get_midpoints(bisection_info.shape, - group, bisection_info, refining_el_old_indices) + midpoints = get_midpoints( + group, tesselation, refining_el_old_indices) new_vertex_indices = np.empty( (new_nelements, group.vertex_indices.shape[1]), @@ -286,17 +142,17 @@ def refine(self, refine_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), + refining_vertices = np.empty(len(tesselation.ref_vertices), dtype=mesh.vertex_id_dtype) refining_vertices.fill(-17) # carry over old vertices - refining_vertices[bisection_info.orig_vertex_indices] = \ + refining_vertices[tesselation.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)): + tesselation.midpoint_indices, + tesselation.midpoint_vertex_pairs)): global_v1 = group.vertex_indices[old_iel, v1] global_v2 = group.vertex_indices[old_iel, v2] @@ -320,7 +176,7 @@ def refine(self, refine_flags): assert (refining_vertices >= 0).all() new_vertex_indices[new_iel_base:new_iel_base+nchildren] = \ - refining_vertices[bisection_info.children] + refining_vertices[tesselation.children] assert (new_vertex_indices >= 0).all() else: @@ -339,8 +195,8 @@ def refine(self, refine_flags): # copy over unchanged nodes new_nodes[:, unrefined_el_new_indices] = group.nodes[:, ~grp_flags] - tesselated_nodes = get_tesselated_nodes(bisection_info.shape, - group, bisection_info, refining_el_old_indices) + tesselated_nodes = get_tesselated_nodes( + group, tesselation, refining_el_old_indices) for old_iel in refining_el_old_indices: new_iel_base = child_el_indices[old_iel] diff --git a/meshmode/mesh/refinement/resampler.py b/meshmode/mesh/refinement/resampler.py deleted file mode 100644 index 39cd539df..000000000 --- a/meshmode/mesh/refinement/resampler.py +++ /dev/null @@ -1,191 +0,0 @@ -__copyright__ = "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. -""" - -from functools import singledispatch - -import numpy as np -import modepy as mp -from pytools import RecordWithoutPickling - -import logging -logger = logging.getLogger(__name__) - - -# {{{ resampling simplex points for refinement - -class SimplexResampler: - @staticmethod - def get_vertex_pair_to_midpoint_order(dim): - return get_vertex_pair_to_midpoint_order(mp.Simplex(dim)) - - @staticmethod - def get_midpoints(group, tesselation, elements): - return get_midpoints(mp.Simplex(group.dim), - group, tesselation, elements) - - @staticmethod - def get_tesselated_nodes(group, tesselation, elements): - return get_tesselated_nodes(mp.Simplex(group.dim), - group, tesselation, elements) - -# }}} - - -# {{{ interface - -# NOTE: internal to refiners: do not make documentation public. - -class TesselationInfo(RecordWithoutPickling): - """ - .. attribute:: ref_vertices - .. attribute:: children - """ - -@singledispatch -def map_unit_nodes_to_children(shape: mp.Shape, tesselation, unit_nodes): - raise NotImplementedError - - -@singledispatch -def get_ref_midpoints(shape: mp.Shape, ref_vertices): - from pytools import add_tuples - space = mp.space_for_shape(shape, 1) - orig_vertices = [ - add_tuples(vt, vt) for vt in mp.node_tuples_for_space(space) - ] - return [rv for rv in ref_vertices if rv not in orig_vertices] - - -@singledispatch -def get_midpoints(shape: mp.Shape, group, tesselation, elements): - """Compute the midpoints of the vertices of the specified elements. - - :arg group: an instance of :class:`meshmode.mesh.MeshElementGroup`. - :arg tesselation: a :class:`TesselationInfo`. - :arg elements: a list of (group-relative) element numbers. - - :return: A :class:`dict` mapping element numbers to midpoint - coordinates, with each value in the map having shape - ``(ambient_dim, nmidpoints)``. The ordering of the midpoints - follows their ordering in the tesselation (see also - :meth:`SimplexResampler.get_vertex_pair_to_midpoint_order`) - """ - if shape.dim != group.dim: - raise ValueError("shape and group dimension do not match") - - if group.vertex_indices is not None: - assert len(group.vertex_indices[0]) == shape.nvertices - - # get midpoints, converted to unit coordinates. - midpoints = -1 + np.array(get_ref_midpoints(shape, tesselation.ref_vertices)) - - space = mp.space_for_shape(shape, group.order) - resampling_mat = mp.resampling_matrix( - mp.basis_for_space(space, shape).functions, - midpoints.T, - group.unit_nodes) - - resampled_midpoints = np.einsum("mu,deu->edm", - resampling_mat, group.nodes[:, elements]) - - return dict(zip(elements, resampled_midpoints)) - - -@singledispatch -def get_tesselated_nodes(shape: mp.Shape, group, tesselation, elements): - """Compute the nodes of the child elements according to the tesselation. - - :arg group: An instance of :class:`meshmode.mesh.MeshElementGroup`. - :arg tesselation: a :class:`TesselationInfo`. - :arg elements: A list of (group-relative) element numbers. - - :return: A :class:`dict` mapping element numbers to node - coordinates, with each value in the map having shape - ``(ambient_dim, nchildren, nunit_nodes)``. - The ordering of the child nodes follows the ordering - of ``tesselation.children.`` - """ - if shape.dim != group.dim: - raise ValueError("shape and group dimension do not match") - - if group.vertex_indices is not None: - assert len(group.vertex_indices[0]) == shape.nvertices - - # get child unit node coordinates. - child_unit_nodes = np.hstack(list( - map_unit_nodes_to_children(shape, tesselation, group.unit_nodes) - )) - - space = mp.space_for_shape(shape, group.order) - resampling_mat = mp.resampling_matrix( - mp.basis_for_space(space, shape).functions, - child_unit_nodes, - group.unit_nodes) - - resampled_unit_nodes = np.einsum("cu,deu->edc", - resampling_mat, group.nodes[:, elements]) - - ambient_dim = len(group.nodes) - nunit_nodes = len(group.unit_nodes[0]) - - return { - el: resampled_unit_nodes[iel].reshape((ambient_dim, -1, nunit_nodes)) - for iel, el in enumerate(elements) - } - raise NotImplementedError(type(shape).__name__) - -# }}} - - -# {{{ simplex - -@map_unit_nodes_to_children.register(mp.Simplex) -def _(shape: mp.Simplex, tesselation, unit_nodes): - ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float) - assert len(unit_nodes.shape) == 2 - - for child_element in tesselation.children: - origin = np.vstack(ref_vertices[child_element[0]]) - basis = ref_vertices.T[:, child_element[1:]] - origin - - yield basis.dot((unit_nodes + 1) / 2) + origin - 1 - -# }}} - - -# {{{ hypercube - -@map_unit_nodes_to_children.register(mp.Hypercube) -def _(shape: mp.Hypercube, tesselation, unit_nodes): - ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float) - assert len(unit_nodes.shape) == 2 - - basis_indices = [1, 2, 4][:shape.dim] - for child_element in tesselation.children: - origin = np.vstack(ref_vertices[child_element[0]]) - basis = ref_vertices.T[:, np.array(child_element)[indices]] - origin - - yield basis.dot((unit_nodes + 1) / 2) + origin - 1 - -# }}} - -# vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index b099fd261..e6868c7d7 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -1,6 +1,7 @@ __copyright__ = """ Copyright (C) 2018 Andreas Kloeckner Copyright (C) 2014-6 Shivam Gupta +Copyright (C) 2020 Alexandru Fikl """ __license__ = """ @@ -23,109 +24,239 @@ THE SOFTWARE. """ +from functools import singledispatch +import numpy as np -from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ - as gnitstam +import modepy as mp +from pytools import RecordWithoutPickling, memoize +import logging +logger = logging.getLogger(__name__) -def mul_tuples(t, m): - return tuple(m * a for a in t) +class TesselationInfo(RecordWithoutPickling): + """ + .. attribute:: children + .. attribute:: ref_vertices -def add_tuples(a, b): - return tuple(ac+bc for ac, bc in zip(a, b)) + .. attribute:: orig_vertex_indices + .. attribute:: midpoint_indices + .. attribute:: midpoint_vertex_pairs + """ -def halve_tuple(a): - def halve(x): - d, r = divmod(x, 2) +class GroupRefinementRecord(RecordWithoutPickling): + """ + .. attribute:: tesselation + .. attribute:: element_mapping + """ + + +def midpoint_tuples(a, b): + def midpoint(x, y): + d, r = divmod(x + y, 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 - - -def tesselatetri(): - result = [] - - node_tuples = list(gnitstam(2, 2)) - node_dict = { - ituple: idx - for idx, ituple in enumerate(node_tuples)} - - def try_add_tri(current, d1, d2, d3): - try: - result.append(( - node_dict[add_tuples(current, d1)], - node_dict[add_tuples(current, d2)], - node_dict[add_tuples(current, d3)], - )) - except KeyError: - pass - - if len(result) > 0: - return [node_tuples, result] - for current in node_tuples: - # this is a tesselation of a square into two triangles. - # subtriangles that fall outside of the master tet are simply not added. - - # 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 - - -def tesselatetet(): - node_tuples = list(gnitstam(2, 3)) - - node_dict = { - ituple: idx - for idx, ituple in enumerate(node_tuples)} - - def try_add_tet(current, d1, d2, d3, d4): - try: - result.append(( - node_dict[add_tuples(current, d1)], - node_dict[add_tuples(current, d2)], - node_dict[add_tuples(current, d3)], - node_dict[add_tuples(current, d4)], - )) - except KeyError: - pass - - result = [] - - if len(result) > 0: - return [node_tuples, result] - for current in node_tuples: - # this is a tesselation of a cube into six tets. - # subtets that fall outside of the master tet are simply not added. - - # positively oriented - try_add_tet(current, (0, 0, 0), (1, 0, 0), (0, 1, 0), (0, 0, 1)) - try_add_tet(current, (1, 0, 1), (1, 0, 0), (0, 0, 1), (0, 1, 0)) - try_add_tet(current, (1, 0, 1), (0, 1, 1), (0, 1, 0), (0, 0, 1)) - - try_add_tet(current, (1, 0, 0), (0, 1, 0), (1, 0, 1), (1, 1, 0)) - 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 - - -def tesselate_simplex_bisection(dim): - if dim == 1: - return tesselateseg() - elif dim == 2: - return tesselatetri() - elif dim == 3: - return tesselatetet() + return tuple(midpoint(ai, bi) for ai, bi in zip(a, b)) + + +# {{{ shape-dependent tesselation helpers + +@singledispatch +def get_child_basis_vertex_indices(shape: mp.Shape, child): + assert len(child) == shape.nvertices + return child[:shape.dim + 1] + + +@get_child_basis_vertex_indices.register(mp.Hypercube) +def _(shape: mp.Hypercube, child): + assert len(child) == shape.nvertices + + # for a cube, we reorder nodes such that the vectors + # (0, 1), (0, 2) and (0, 4) + # form an orthogonal basis, since (0, 3) is linearly dependent on 1 and 2. + + return [child[i] for i in [0, 1, 2, 4][:shape.dim + 1]] + +# }}} + + +# {{{ resampling + +def get_ref_midpoints(shape, ref_vertices): + r"""The reference element is considered to be, e.g. for a 2 simplex:: + + F + | \ + | \ + D----E + | /| \ + | / | \ + A----B----C + + where the midpoints are ``(B, E, D)``. The same applies to other shapes + and higher dimensions. + + :arg ref_vertices: a :class:`list` of node index :class:`tuple`\ s + on :math:`[0, 2]^d`. + """ + + from pytools import add_tuples + space = mp.space_for_shape(shape, 1) + orig_vertices = [ + add_tuples(vt, vt) for vt in mp.node_tuples_for_space(space) + ] + return [rv for rv in ref_vertices if rv not in orig_vertices] + + +def map_unit_nodes_to_children(shape, unit_nodes, tesselation): + """ + :arg unit_nodes: an :class:`~numpy.ndarray` of shape ``(dim, nnodes)``. + :arg tesselation: a :class:`TesselationInfo`. + """ + + ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float).T + assert len(unit_nodes.shape) == 2 + + for child_element in tesselation.children: + indices = get_child_basis_vertex_indices(shape, child_element) + + origin = ref_vertices[:, indices[0]].reshape(-1, 1) + basis = ref_vertices[:, indices[1:]] - origin + + # mapped nodes are on [0, 2], so we subtract 1 to get it to [-1, 1] + yield basis.dot((unit_nodes + 1.0) / 2.0) + origin - 1.0 + + +def get_midpoints(group, tesselation, elements): + """Compute the midpoints of the vertices of the specified elements. + + :arg group: an instance of :class:`meshmode.mesh.MeshElementGroup`. + :arg tesselation: a :class:`TesselationInfo`. + :arg elements: a list of (group-relative) element numbers. + + :return: A :class:`dict` mapping element numbers to midpoint + coordinates, with each value in the map having shape + ``(ambient_dim, nmidpoints)``. The ordering of the midpoints + follows their ordering in the tesselation. + """ + from meshmode.mesh import _ModepyElementGroup + if not isinstance(group, _ModepyElementGroup): + raise TypeError(f"groups of type '{type(group.mesh_el_group).__name__}'" + " are not supported.") + + shape = group._modepy_shape + space = mp.space_for_shape(shape, group.order) + + # get midpoints in reference coordinates + midpoints = -1 + np.array(get_ref_midpoints(shape, tesselation.ref_vertices)) + + # resample midpoints to ambient coordinates + resampling_mat = mp.resampling_matrix( + mp.basis_for_space(space, shape).functions, + midpoints.T, + group.unit_nodes) + + resampled_midpoints = np.einsum("mu,deu->edm", + resampling_mat, group.nodes[:, elements]) + + return dict(zip(elements, resampled_midpoints)) + + +def get_tesselated_nodes(group, tesselation, elements): + """Compute the nodes of the child elements according to the tesselation. + + :arg group: An instance of :class:`meshmode.mesh.MeshElementGroup`. + :arg tesselation: a :class:`TesselationInfo`. + :arg elements: A list of (group-relative) element numbers. + + :return: A :class:`dict` mapping element numbers to node + coordinates, with each value in the map having shape + ``(ambient_dim, nchildren, nunit_nodes)``. + The ordering of the child nodes follows the ordering + of ``tesselation.children.`` + """ + from meshmode.mesh import _ModepyElementGroup + if not isinstance(group, _ModepyElementGroup): + raise TypeError(f"groups of type '{type(group.mesh_el_group).__name__}'" + " are not supported.") + + shape = group._modepy_shape + space = mp.space_for_shape(shape, group.order) + + # get child unit node coordinates. + child_unit_nodes = np.hstack(list( + map_unit_nodes_to_children(shape, group.unit_nodes, tesselation) + )) + + # resample child nodes to ambient coordinates + resampling_mat = mp.resampling_matrix( + mp.basis_for_space(space, shape).functions, + child_unit_nodes, + group.unit_nodes) + + resampled_unit_nodes = np.einsum("cu,deu->edc", + resampling_mat, group.nodes[:, elements]) + + ambient_dim = len(group.nodes) + nunit_nodes = len(group.unit_nodes[0]) + + return { + el: resampled_unit_nodes[iel].reshape((ambient_dim, -1, nunit_nodes)) + for iel, el in enumerate(elements) + } + +# }}} + + +# {{{ tesselation + +def get_tesselation_info(shape): + space = mp.space_for_shape(shape, 2) + ref_vertices = mp.node_tuples_for_space(space) + ref_vertices_to_index = {rv: i for i, rv in enumerate(ref_vertices)} + + from pytools import add_tuples + space = mp.space_for_shape(shape, 1) + orig_vertices = tuple([ + add_tuples(vt, vt) for vt in mp.node_tuples_for_space(space) + ]) + orig_vertex_indices = [ref_vertices_to_index[vt] for vt in orig_vertices] + + midpoints = get_ref_midpoints(shape, ref_vertices) + midpoint_indices = [ref_vertices_to_index[mp] for mp in midpoints] + + midpoint_to_vertex_pairs = { + midpoint: (i, j) + for i, ivt in enumerate(orig_vertices) + for j, jvt in enumerate(orig_vertices) + for midpoint in [midpoint_tuples(ivt, jvt)] + if i < j and midpoint in midpoints + } + # ensure order matches the one in midpoint_indices + midpoint_vertex_pairs = [midpoint_to_vertex_pairs[m] for m in midpoints] + + return TesselationInfo( + ref_vertices=ref_vertices, + children=np.array(mp.submesh_for_shape(shape, ref_vertices)), + orig_vertex_indices=np.array(orig_vertex_indices), + midpoint_indices=np.array(midpoint_indices), + midpoint_vertex_pairs=midpoint_vertex_pairs, + ) + + +@memoize +def get_bisection_tesselation_info(group_type, dim): + from meshmode.mesh import _ModepyElementGroup + if issubclass(group_type, _ModepyElementGroup): + from meshmode.mesh.refinement.tesselate import get_tesselation_info + shape = group_type._modepy_shape_cls(dim) + return get_tesselation_info(shape) else: - raise ValueError("cannot tesselate %d-simplex" % dim) + raise NotImplementedError( + "bisection for element groups of type {group_type.__name__}") + +# }}} diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py deleted file mode 100644 index 1b717a803..000000000 --- a/meshmode/mesh/refinement/utils.py +++ /dev/null @@ -1,160 +0,0 @@ -__copyright__ = "Copyright (C) 2014-6 Shivam Gupta, Andreas Kloeckner" - -__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 - -import logging -logger = logging.getLogger(__name__) - - -# {{{ map unit nodes to children - - -def map_unit_nodes_to_children(unit_nodes, tesselation): - """ - Given a collection of unit nodes, return the coordinates of the - unit nodes mapped onto each of the children of the reference - element. - - The tesselation should follow the format of - :func:`meshmode.mesh.tesselate.tesselatetri()` or - :func:`meshmode.mesh.tesselate.tesselatetet()`. - - `unit_nodes` should be relative to the unit simplex coordinates in - :module:`modepy`. - - :arg unit_nodes: shaped `(dim, nunit_nodes)` - :arg tesselation: With attributes `ref_vertices`, `children` - """ - ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float) - assert len(unit_nodes.shape) == 2 - - for child_element in tesselation.children: - center = np.vstack(ref_vertices[child_element[0]]) - # Scale by 1/2 since sides in the tesselation have length 2. - aff_mat = (ref_vertices.T[:, child_element[1:]] - center) / 2 - # (-1, -1, ...) in unit_nodes = (0, 0, ...) in ref_vertices. - # Hence the translation by +/- 1. - yield aff_mat.dot(unit_nodes + 1) + center - 1 - -# }}} - -# {{{ test nodal adjacency against geometry - - -def is_symmetric(relation, debug=False): - for a, other_list in enumerate(relation): - for b in other_list: - if a not in relation[b]: - if debug: - print("Relation is not symmetric: %s -> %s, but not %s -> %s" - % (a, b, b, a)) - return False - - return True - - -def check_nodal_adj_against_geometry(mesh, tol=1e-12): - def group_and_iel_to_global_iel(igrp, iel): - return mesh.groups[igrp].element_nr_base + iel - - logger.debug("nodal adj test: tree build") - from meshmode.mesh.tools import make_element_lookup_tree - tree = make_element_lookup_tree(mesh, eps=tol) - logger.debug("nodal adj test: tree build done") - - from meshmode.mesh.processing import find_bounding_box - bbox_min, bbox_max = find_bounding_box(mesh) - - nadj = mesh.nodal_adjacency - nvertices_per_element = len(mesh.groups[0].vertex_indices[0]) - - connected_to_element_geometry = [set() for i in range(mesh.nelements)] - connected_to_element_connectivity = [set() for i in range(mesh.nelements)] - - for igrp, grp in enumerate(mesh.groups): - for iel_grp in range(grp.nelements): - iel_g = group_and_iel_to_global_iel(igrp, iel_grp) - nb_starts = nadj.neighbors_starts - for nb_iel_g in nadj.neighbors[nb_starts[iel_g]:nb_starts[iel_g+1]]: - connected_to_element_connectivity[iel_g].add(nb_iel_g) - - for vertex_index in grp.vertex_indices[iel_grp]: - vertex = mesh.vertices[:, vertex_index] - - # check which elements touch this vertex - for nearby_igrp, nearby_iel in tree.generate_matches(vertex): - if nearby_igrp == igrp and nearby_iel == iel_grp: - continue - nearby_grp = mesh.groups[nearby_igrp] - - nearby_origin_vertex = mesh.vertices[ - :, nearby_grp.vertex_indices[nearby_iel][0]] # noqa - transformation = np.empty( - (len(mesh.vertices), nvertices_per_element-1)) - vertex_transformed = vertex - nearby_origin_vertex - - for inearby_vertex_index, nearby_vertex_index in enumerate( - nearby_grp.vertex_indices[nearby_iel][1:]): - nearby_vertex = mesh.vertices[:, nearby_vertex_index] - transformation[:, inearby_vertex_index] = \ - nearby_vertex - nearby_origin_vertex - bary_coord, residual = \ - np.linalg.lstsq(transformation, vertex_transformed)[0:2] - - is_in_element_span = ( - residual.size == 0 - or np.linalg.norm(vertex_transformed) == 0 - or (np.linalg.norm(residual) - / np.linalg.norm(vertex_transformed)) <= tol) - - is_connected = ( - is_in_element_span - and np.sum(bary_coord) <= 1+tol - and (bary_coord >= -tol).all()) - el1 = group_and_iel_to_global_iel(nearby_igrp, nearby_iel) - el2 = group_and_iel_to_global_iel(igrp, iel_grp) - - if is_connected: - connected_to_element_geometry[el1].add(el2) - connected_to_element_geometry[el2].add(el1) - - assert is_symmetric(connected_to_element_connectivity, debug=True) - - # The geometric adjacency relation isn't necessary symmetric: - # - # /| - # / | - # / |\ - # B \ |/ A - # \ | - # \| - # - # Element A will see element B (its vertices are near B) but not the other - # way around. - - assert connected_to_element_geometry == connected_to_element_connectivity - -# }}} - -# vim: foldmethod=marker From 33ddc28c7c8f45dfff0cf582872cb2fd3eab5b3e Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 6 Dec 2020 11:21:11 -0600 Subject: [PATCH 03/20] add test for tensor product refinement connection --- meshmode/mesh/refinement/__init__.py | 1 + meshmode/mesh/refinement/tesselate.py | 10 +- meshmode/mesh/refinement/utils.py | 127 ++++++++++++++++++++++++++ test/test_refinement.py | 40 +++++--- 4 files changed, 161 insertions(+), 17 deletions(-) create mode 100644 meshmode/mesh/refinement/utils.py diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index 53ecbd162..bb56f21e2 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -625,6 +625,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr nelements_in_grp = grp.nelements del self.group_refinement_records[:] + from meshmode.mesh import SimplexElementGroup for grp_idx, grp in enumerate(self.last_mesh.groups): if not isinstance(grp, SimplexElementGroup): raise TypeError("refinement not supported for groups of type " diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index e6868c7d7..08fdf344a 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -75,9 +75,13 @@ def get_child_basis_vertex_indices(shape: mp.Shape, child): def _(shape: mp.Hypercube, child): assert len(child) == shape.nvertices - # for a cube, we reorder nodes such that the vectors - # (0, 1), (0, 2) and (0, 4) - # form an orthogonal basis, since (0, 3) is linearly dependent on 1 and 2. + # * for a cube, we reorder nodes such that the vectors + # + # (0, 1), (0, 2) and (0, 4) + # + # form an orthogonal basis, since (0, 3) is linearly dependent on 1 and 2. + # + # * lines and squares don't require any reordering return [child[i] for i in [0, 1, 2, 4][:shape.dim + 1]] diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py new file mode 100644 index 000000000..446dbea02 --- /dev/null +++ b/meshmode/mesh/refinement/utils.py @@ -0,0 +1,127 @@ +__copyright__ = "Copyright (C) 2014-6 Shivam Gupta, Andreas Kloeckner" + +__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 + +import logging +logger = logging.getLogger(__name__) + + +# {{{ test nodal adjacency against geometry + +def is_symmetric(relation, debug=False): + for a, other_list in enumerate(relation): + for b in other_list: + if a not in relation[b]: + if debug: + print("Relation is not symmetric: %s -> %s, but not %s -> %s" + % (a, b, b, a)) + return False + + return True + + +def check_nodal_adj_against_geometry(mesh, tol=1e-12): + def group_and_iel_to_global_iel(igrp, iel): + return mesh.groups[igrp].element_nr_base + iel + + logger.debug("nodal adj test: tree build") + from meshmode.mesh.tools import make_element_lookup_tree + tree = make_element_lookup_tree(mesh, eps=tol) + logger.debug("nodal adj test: tree build done") + + from meshmode.mesh.processing import find_bounding_box + bbox_min, bbox_max = find_bounding_box(mesh) + + nadj = mesh.nodal_adjacency + nvertices_per_element = len(mesh.groups[0].vertex_indices[0]) + + connected_to_element_geometry = [set() for i in range(mesh.nelements)] + connected_to_element_connectivity = [set() for i in range(mesh.nelements)] + + for igrp, grp in enumerate(mesh.groups): + for iel_grp in range(grp.nelements): + iel_g = group_and_iel_to_global_iel(igrp, iel_grp) + nb_starts = nadj.neighbors_starts + for nb_iel_g in nadj.neighbors[nb_starts[iel_g]:nb_starts[iel_g+1]]: + connected_to_element_connectivity[iel_g].add(nb_iel_g) + + for vertex_index in grp.vertex_indices[iel_grp]: + vertex = mesh.vertices[:, vertex_index] + + # check which elements touch this vertex + for nearby_igrp, nearby_iel in tree.generate_matches(vertex): + if nearby_igrp == igrp and nearby_iel == iel_grp: + continue + nearby_grp = mesh.groups[nearby_igrp] + + nearby_origin_vertex = mesh.vertices[ + :, nearby_grp.vertex_indices[nearby_iel][0]] # noqa + transformation = np.empty( + (len(mesh.vertices), nvertices_per_element-1)) + vertex_transformed = vertex - nearby_origin_vertex + + for inearby_vertex_index, nearby_vertex_index in enumerate( + nearby_grp.vertex_indices[nearby_iel][1:]): + nearby_vertex = mesh.vertices[:, nearby_vertex_index] + transformation[:, inearby_vertex_index] = \ + nearby_vertex - nearby_origin_vertex + bary_coord, residual = \ + np.linalg.lstsq(transformation, vertex_transformed)[0:2] + + is_in_element_span = ( + residual.size == 0 + or np.linalg.norm(vertex_transformed) == 0 + or (np.linalg.norm(residual) + / np.linalg.norm(vertex_transformed)) <= tol) + + is_connected = ( + is_in_element_span + and np.sum(bary_coord) <= 1+tol + and (bary_coord >= -tol).all()) + el1 = group_and_iel_to_global_iel(nearby_igrp, nearby_iel) + el2 = group_and_iel_to_global_iel(igrp, iel_grp) + + if is_connected: + connected_to_element_geometry[el1].add(el2) + connected_to_element_geometry[el2].add(el1) + + assert is_symmetric(connected_to_element_connectivity, debug=True) + + # The geometric adjacency relation isn't necessary symmetric: + # + # /| + # / | + # / |\ + # B \ |/ A + # \ | + # \| + # + # Element A will see element B (its vertices are near B) but not the other + # way around. + + assert connected_to_element_geometry == connected_to_element_connectivity + +# }}} + +# vim: foldmethod=marker diff --git a/test/test_refinement.py b/test/test_refinement.py index 9314c93fa..c88f7a148 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -26,12 +26,13 @@ import numpy as np import pytest -from meshmode.array_context import ( # noqa +from meshmode import _acf # noqa: F401 +from meshmode.array_context import ( # noqa: F401 pytest_generate_tests_for_pyopencl_array_context as pytest_generate_tests) from meshmode.dof_array import thaw -from meshmode.mesh.generation import ( # noqa +from meshmode.mesh.generation import ( # noqa: F401 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, RefinerWithoutAdjacency @@ -41,6 +42,8 @@ InterpolatoryQuadratureSimplexGroupFactory, PolynomialWarpAndBlendGroupFactory, PolynomialEquidistantSimplexGroupFactory, + LegendreGaussLobattoTensorProductGroupFactory, + GaussLegendreTensorProductGroupFactory, ) logger = logging.getLogger(__name__) @@ -156,14 +159,17 @@ 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, - PolynomialEquidistantSimplexGroupFactory +@pytest.mark.parametrize(("refiner_cls", "group_factory"), [ + (Refiner, InterpolatoryQuadratureSimplexGroupFactory), + (Refiner, PolynomialWarpAndBlendGroupFactory), + (Refiner, PolynomialEquidistantSimplexGroupFactory), + + (RefinerWithoutAdjacency, InterpolatoryQuadratureSimplexGroupFactory), + (RefinerWithoutAdjacency, PolynomialWarpAndBlendGroupFactory), + (RefinerWithoutAdjacency, PolynomialEquidistantSimplexGroupFactory), + + (RefinerWithoutAdjacency, LegendreGaussLobattoTensorProductGroupFactory), + (RefinerWithoutAdjacency, GaussLegendreTensorProductGroupFactory), ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ ("circle", 1, [20, 30, 40]), @@ -182,14 +188,19 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): def test_refinement_connection( actx_factory, refiner_cls, group_factory, mesh_name, dim, mesh_pars, mesh_order, refine_flags, visualize=False): + group_cls = group_factory.mesh_group_class + if issubclass(group_cls, TensorProductElementGroup): + if mesh_name in ["circle", "blob"]: + pytest.skip("mesh does not have tensor product support") + from random import seed seed(13) - # Discretization order - order = 5 - actx = actx_factory() + # discretization order + order = 5 + from meshmode.discretization import Discretization from meshmode.discretization.connection import ( make_refinement_connection, check_connection) @@ -214,7 +225,8 @@ def test_refinement_connection( h = float(mesh_par) elif mesh_name == "warp": from meshmode.mesh.generation import generate_warped_rect_mesh - mesh = generate_warped_rect_mesh(dim, order=mesh_order, n=mesh_par) + mesh = generate_warped_rect_mesh(dim, order=mesh_order, n=mesh_par, + group_cls=group_cls) h = 1/mesh_par else: raise ValueError("mesh_name not recognized") From 85d8762a5c37c6f2e7411203de52fbe8e20d3c3a Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Dec 2020 21:29:08 -0600 Subject: [PATCH 04/20] use dataclasses and add docs --- meshmode/mesh/refinement/__init__.py | 1 - meshmode/mesh/refinement/no_adjacency.py | 3 +- meshmode/mesh/refinement/tesselate.py | 49 ++++++++++++++++++++++-- 3 files changed, 46 insertions(+), 7 deletions(-) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index bb56f21e2..278ba08fe 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -23,7 +23,6 @@ import numpy as np import itertools -from pytools import RecordWithoutPickling from meshmode.mesh.refinement.no_adjacency import ( # noqa: F401 RefinerWithoutAdjacency) diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index 484363c72..569c2e298 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -26,8 +26,6 @@ import numpy as np -from meshmode.mesh.refinement.tesselate import GroupRefinementRecord - import logging logger = logging.getLogger(__name__) @@ -115,6 +113,7 @@ def refine(self, refine_flags): # }}} + from meshmode.mesh.refinement.tesselate import GroupRefinementRecord group_refinement_records.append( GroupRefinementRecord( tesselation=tesselation, diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index 08fdf344a..81e7c48e6 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -24,33 +24,74 @@ THE SOFTWARE. """ +from dataclasses import dataclass from functools import singledispatch -import numpy as np +import numpy as np import modepy as mp -from pytools import RecordWithoutPickling, memoize + +from pytools import memoize import logging logger = logging.getLogger(__name__) +from typing import List, Tuple, Optional -class TesselationInfo(RecordWithoutPickling): + +@dataclass(frozen=True) +class TesselationInfo: """ .. attribute:: children + + A tesselation of the reference element, given here by + :attr:`ref_vertices`. + .. attribute:: ref_vertices + A list of tuples (similar to :func:`modepy.node_tuples_for_space`) + for the reference element containing midpoints. This is equivalent + to a second-order equidistant element. + .. attribute:: orig_vertex_indices + + Indices into :attr:`ref_vertices` that select only the vertices, i.e. + without the midpoints. + .. attribute:: midpoint_indices + + Indices into :attr:`ref_vertices` that select only the midpoints, i.e. + without :attr:`orig_vertex_indices`. + .. attribute:: midpoint_vertex_pairs + + A list of tuples ``(v1, v2)`` of indices into :attr:`orig_vertex_indices` + that give for each midpoint the two vertices on the same line. """ + children: np.ndarray + ref_vertices: List[Tuple[int, ...]] + + orig_vertex_indices: Optional[np.ndarray] = None + midpoint_indices: Optional[np.ndarray] = None + midpoint_vertex_pairs: Optional[List[Tuple[int, int]]] = None + -class GroupRefinementRecord(RecordWithoutPickling): +@dataclass(frozen=True) +class GroupRefinementRecord: """ .. attribute:: tesselation + + A :class:`TesselationInfo` that describes the tesselation of the + element group. + .. attribute:: element_mapping + + A mapping from the original elements to the refined child elements. """ + tesselation: TesselationInfo + element_mapping: List[List[int]] + def midpoint_tuples(a, b): def midpoint(x, y): From 13de25e7d2132057d373e425158769a2146fc7a1 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Dec 2020 21:35:25 -0600 Subject: [PATCH 05/20] remove unused functions --- meshmode/mesh/refinement/__init__.py | 30 +++++++++++----------------- 1 file changed, 12 insertions(+), 18 deletions(-) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index 278ba08fe..4716fdbba 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -20,12 +20,12 @@ THE SOFTWARE. """ +import itertools +from functools import partial import numpy as np -import itertools -from meshmode.mesh.refinement.no_adjacency import ( # noqa: F401 - RefinerWithoutAdjacency) +from meshmode.mesh.refinement.no_adjacency import RefinerWithoutAdjacency from meshmode.mesh.refinement.tesselate import \ TesselationInfo, GroupRefinementRecord @@ -38,6 +38,11 @@ .. autofunction :: refine_uniformly """ +__all__ = [ + "Refiner", "RefinerWithoutAdjacency", "refine_uniformly" +] + + # {{{ deprecated class SimplexResampler: @@ -68,20 +73,9 @@ def tesselate_simplex(dim): node_tuples = mp.node_tuples_for_space(space) return node_tuples, mp.submesh_for_shape(shape, node_tuples) - -def tesselateseg(): - return tesselate_simplex(1) - - -def tesselatetri(): - return tesselate_simplex(2) - - -def tesselatetet(): - return tesselate_simplex(3) - # }}} + class TreeRayNode: """Describes a ray as a tree, this class represents each node in this tree .. attribute:: left @@ -138,9 +132,9 @@ def __init__(self, mesh): self.lazy = False self.seen_tuple = {} self.group_refinement_records = [] - seg_node_tuples, seg_result = tesselateseg() - tri_node_tuples, tri_result = tesselatetri() - tet_node_tuples, tet_result = tesselatetet() + seg_node_tuples, seg_result = tesselate_simplex(1) + tri_node_tuples, tri_result = tesselate_simplex(2) + tet_node_tuples, tet_result = tesselate_simplex(3) #quadrilateral_node_tuples = [ #print tri_result, tet_result self.simplex_node_tuples = [ From ac1e905bea4620ed2cfe816aaff59be1b33be6a5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 11 Dec 2020 21:41:28 -0600 Subject: [PATCH 06/20] remove unused import and improve naming --- meshmode/mesh/refinement/__init__.py | 4 ++-- meshmode/mesh/refinement/no_adjacency.py | 12 ++++++------ meshmode/mesh/refinement/tesselate.py | 11 +++++------ 3 files changed, 13 insertions(+), 14 deletions(-) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index 4716fdbba..08c1c6b0e 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -57,12 +57,12 @@ def get_vertex_pair_to_midpoint_order(dim): @staticmethod def get_midpoints(group, tesselation, elements): import meshmode.mesh.refinement.tesselate as tess - return tess.get_midpoints(group, tesselation, elements) + return tess.get_group_midpoints(group, tesselation, elements) @staticmethod def get_tesselated_nodes(group, tesselation, elements): import meshmode.mesh.refinement.tesselate as tess - return tess.get_tesselated_nodes(group, tesselation, elements) + return tess.get_group_tesselated_nodes(group, tesselation, elements) def tesselate_simplex(dim): diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index 569c2e298..0376dcdb8 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -84,12 +84,12 @@ def refine(self, refine_flags): inew_vertex = mesh.nvertices from meshmode.mesh.refinement.tesselate import ( - get_bisection_tesselation_info, - get_midpoints, - get_tesselated_nodes) + get_group_tesselation_info, + get_group_midpoints, + get_group_tesselated_nodes) for igrp, group in enumerate(mesh.groups): - tesselation = get_bisection_tesselation_info(type(group), group.dim) + tesselation = get_group_tesselation_info(type(group), group.dim) # {{{ compute counts and index arrays @@ -126,7 +126,7 @@ def refine(self, refine_flags): # {{{ get new vertices together if perform_vertex_updates: - midpoints = get_midpoints( + midpoints = get_group_midpoints( group, tesselation, refining_el_old_indices) new_vertex_indices = np.empty( @@ -194,7 +194,7 @@ def refine(self, refine_flags): # copy over unchanged nodes new_nodes[:, unrefined_el_new_indices] = group.nodes[:, ~grp_flags] - tesselated_nodes = get_tesselated_nodes( + tesselated_nodes = get_group_tesselated_nodes( group, tesselation, refining_el_old_indices) for old_iel in refining_el_old_indices: diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index 81e7c48e6..1a57e75d6 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -176,7 +176,7 @@ def map_unit_nodes_to_children(shape, unit_nodes, tesselation): yield basis.dot((unit_nodes + 1.0) / 2.0) + origin - 1.0 -def get_midpoints(group, tesselation, elements): +def get_group_midpoints(group, tesselation, elements): """Compute the midpoints of the vertices of the specified elements. :arg group: an instance of :class:`meshmode.mesh.MeshElementGroup`. @@ -211,7 +211,7 @@ def get_midpoints(group, tesselation, elements): return dict(zip(elements, resampled_midpoints)) -def get_tesselated_nodes(group, tesselation, elements): +def get_group_tesselated_nodes(group, tesselation, elements): """Compute the nodes of the child elements according to the tesselation. :arg group: An instance of :class:`meshmode.mesh.MeshElementGroup`. @@ -259,7 +259,7 @@ def get_tesselated_nodes(group, tesselation, elements): # {{{ tesselation -def get_tesselation_info(shape): +def get_shape_tesselation_info(shape): space = mp.space_for_shape(shape, 2) ref_vertices = mp.node_tuples_for_space(space) ref_vertices_to_index = {rv: i for i, rv in enumerate(ref_vertices)} @@ -294,12 +294,11 @@ def get_tesselation_info(shape): @memoize -def get_bisection_tesselation_info(group_type, dim): +def get_group_tesselation_info(group_type, dim): from meshmode.mesh import _ModepyElementGroup if issubclass(group_type, _ModepyElementGroup): - from meshmode.mesh.refinement.tesselate import get_tesselation_info shape = group_type._modepy_shape_cls(dim) - return get_tesselation_info(shape) + return get_shape_tesselation_info(shape) else: raise NotImplementedError( "bisection for element groups of type {group_type.__name__}") From c3d7335529e95bf6179fc385c33e671da5f20fb7 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Wed, 16 Dec 2020 13:57:38 -0600 Subject: [PATCH 07/20] add vim foldmethod MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- meshmode/mesh/refinement/tesselate.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index 1a57e75d6..54c980387 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -304,3 +304,5 @@ def get_group_tesselation_info(group_type, dim): "bisection for element groups of type {group_type.__name__}") # }}} + +# vim: foldmethod=marker From b7e1f8f102982e7950299f2c63f8a8183b41a37d Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 16 Dec 2020 14:06:41 -0600 Subject: [PATCH 08/20] deprecate Refiner --- meshmode/mesh/refinement/__init__.py | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index 08c1c6b0e..a9ed8050a 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -123,6 +123,10 @@ class Refiner: # {{{ constructor def __init__(self, mesh): + from warnings import warn + warn("Refiner is deprecated and will be removed in 2022.", + DeprecationWarning, stacklevel=2) + 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 " @@ -974,9 +978,13 @@ def generate_nodal_adjacency(self, nelements, nvertices, groups): def refine_uniformly(mesh, iterations, with_adjacency=False): if with_adjacency: - refiner = Refiner(mesh) - else: - refiner = RefinerWithoutAdjacency(mesh) + # For conforming meshes, even RefinerWithoutAdjacency will reconstruct + # adjacency from vertex identity. + + if not mesh.is_conforming: + raise ValueError("mesh must be conforming if adjacency is desired") + + refiner = RefinerWithoutAdjacency(mesh) for _ in range(iterations): refiner.refine_uniformly() From e3cb2970b0c0d8ac4b9d15bfdb620d5525b85d9f Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 16 Dec 2020 14:51:08 -0600 Subject: [PATCH 09/20] fix and document map_unit_nodes_to_children --- .../discretization/connection/refinement.py | 11 +--- meshmode/mesh/refinement/tesselate.py | 47 +-------------- meshmode/mesh/refinement/utils.py | 58 +++++++++++++++++++ 3 files changed, 63 insertions(+), 53 deletions(-) diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index f09683050..df31c8f78 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -89,14 +89,9 @@ def _build_interpolation_batches_for_group( fine_unit_nodes = fine_discr_group.unit_nodes meg = fine_discr_group.mesh_el_group - from meshmode.mesh import _ModepyElementGroup - if isinstance(meg, _ModepyElementGroup): - from meshmode.mesh.refinement.tesselate import map_unit_nodes_to_children - mapped_unit_nodes = map_unit_nodes_to_children( - meg._modepy_shape, fine_unit_nodes, record.tesselation) - else: - raise TypeError("cannot construct refinement connection for groups " - f"of type '{type(meg).__name__}'") + from meshmode.mesh.refinement.utils import map_unit_nodes_to_children + mapped_unit_nodes = map_unit_nodes_to_children( + meg, fine_unit_nodes, record.tesselation) from itertools import chain for from_bin, to_bin, unit_nodes in zip( diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index 54c980387..d1fad5871 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -104,31 +104,6 @@ def midpoint(x, y): return tuple(midpoint(ai, bi) for ai, bi in zip(a, b)) -# {{{ shape-dependent tesselation helpers - -@singledispatch -def get_child_basis_vertex_indices(shape: mp.Shape, child): - assert len(child) == shape.nvertices - return child[:shape.dim + 1] - - -@get_child_basis_vertex_indices.register(mp.Hypercube) -def _(shape: mp.Hypercube, child): - assert len(child) == shape.nvertices - - # * for a cube, we reorder nodes such that the vectors - # - # (0, 1), (0, 2) and (0, 4) - # - # form an orthogonal basis, since (0, 3) is linearly dependent on 1 and 2. - # - # * lines and squares don't require any reordering - - return [child[i] for i in [0, 1, 2, 4][:shape.dim + 1]] - -# }}} - - # {{{ resampling def get_ref_midpoints(shape, ref_vertices): @@ -157,25 +132,6 @@ def get_ref_midpoints(shape, ref_vertices): return [rv for rv in ref_vertices if rv not in orig_vertices] -def map_unit_nodes_to_children(shape, unit_nodes, tesselation): - """ - :arg unit_nodes: an :class:`~numpy.ndarray` of shape ``(dim, nnodes)``. - :arg tesselation: a :class:`TesselationInfo`. - """ - - ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float).T - assert len(unit_nodes.shape) == 2 - - for child_element in tesselation.children: - indices = get_child_basis_vertex_indices(shape, child_element) - - origin = ref_vertices[:, indices[0]].reshape(-1, 1) - basis = ref_vertices[:, indices[1:]] - origin - - # mapped nodes are on [0, 2], so we subtract 1 to get it to [-1, 1] - yield basis.dot((unit_nodes + 1.0) / 2.0) + origin - 1.0 - - def get_group_midpoints(group, tesselation, elements): """Compute the midpoints of the vertices of the specified elements. @@ -233,8 +189,9 @@ def get_group_tesselated_nodes(group, tesselation, elements): space = mp.space_for_shape(shape, group.order) # get child unit node coordinates. + from meshmode.mesh.refinement.utils import map_unit_nodes_to_children child_unit_nodes = np.hstack(list( - map_unit_nodes_to_children(shape, group.unit_nodes, tesselation) + map_unit_nodes_to_children(group, group.unit_nodes, tesselation) )) # resample child nodes to ambient coordinates diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index 446dbea02..1e31f1ab8 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -20,13 +20,71 @@ THE SOFTWARE. """ +from functools import singledispatch import numpy as np +from meshmode.mesh import ( + MeshElementGroup, + SimplexElementGroup, + TensorProductElementGroup) + import logging logger = logging.getLogger(__name__) +# {{{ map child unit nodes + +@singledispatch +def map_unit_nodes_to_children(meg: MeshElementGroup, + unit_nodes, tess_info) -> np.ndarray: + """ + :arg unit_nodes: an :class:`~numpy.ndarray` of unit nodes on the + element type described by *meg*. + :arg tess_info: a :class:`~meshmode.mesh.refinement.tesselate.TesselationInfo`. + :returns: an :class:`~numpy.ndarray` of mapped unit nodes for each + child in the tesselation. + """ + raise NotImplementedError(type(meg).__name__) + + +@singledispatch.register(SimplexElementGroup) +def _(meg: SimplexElementGroup, unit_nodes, tess_info): + ref_vertices = np.array(tess_info.ref_vertices, dtype=np.float).T + assert len(unit_nodes.shape) == 2 + + for child in tess_info.children: + origin = ref_vertices[:, child[0]].reshape(-1, 1) + basis = ref_vertices[:, child[1:]] - origin + + # mapped nodes are on [0, 2], so we subtract 1 to get it to [-1, 1] + yield basis.dot((unit_nodes + 1.0) / 2.0) + origin - 1.0 + + +@singledispatch.register(TensorProductElementGroup) +def _(meg: TensorProductElementGroup, unit_nodes, tess_info): + ref_vertices = np.array(tess_info.ref_vertices, dtype=np.float).T + assert len(unit_nodes.shape) == 2 + + import modepy as mp + from dataclass import replace + space = replace(meg._modepy_space, order=1) + basis_indices = [ + i for i, node in mp.node_tuples_for_space(space) + if sum(node) == 1 + ] + + for child in tess_info.children: + child_arr = np.array(child) + origin = ref_vertices[:, child[0]].reshape(-1, 1) + basis = ref_vertices[:, child[basis_indices]] - origin + + # mapped nodes are on [0, 2], so we subtract 1 to get it to [-1, 1] + yield basis.dot((unit_nodes + 1.0) / 2.0) + origin - 1.0 + +# }}} + + # {{{ test nodal adjacency against geometry def is_symmetric(relation, debug=False): From 0faf929b8575446b0f66da6d016f9907f6ec4f9c Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 16 Dec 2020 15:24:46 -0600 Subject: [PATCH 10/20] singledispatch more in refinement --- meshmode/mesh/refinement/no_adjacency.py | 2 +- meshmode/mesh/refinement/tesselate.py | 146 ++++++++++++----------- meshmode/mesh/refinement/utils.py | 17 +-- 3 files changed, 83 insertions(+), 82 deletions(-) diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index 0376dcdb8..dda338b08 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -89,7 +89,7 @@ def refine(self, refine_flags): get_group_tesselated_nodes) for igrp, group in enumerate(mesh.groups): - tesselation = get_group_tesselation_info(type(group), group.dim) + tesselation = get_group_tesselation_info(group) # {{{ compute counts and index arrays diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index d1fad5871..43d85c1df 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -28,9 +28,9 @@ from functools import singledispatch import numpy as np -import modepy as mp -from pytools import memoize +import modepy as mp +from meshmode.mesh import MeshElementGroup, _ModepyElementGroup import logging logger = logging.getLogger(__name__) @@ -38,6 +38,8 @@ from typing import List, Tuple, Optional +# {{{ interface + @dataclass(frozen=True) class TesselationInfo: """ @@ -93,7 +95,52 @@ class GroupRefinementRecord: element_mapping: List[List[int]] -def midpoint_tuples(a, b): +@singledispatch +def get_group_midpoints(meg: MeshElementGroup, tesselation, elements): + """Compute the midpoints of the vertices of the specified elements. + + :arg group: an instance of :class:`meshmode.mesh.MeshElementGroup`. + :arg tesselation: a :class:`TesselationInfo`. + :arg elements: a list of (group-relative) element numbers. + + :return: A :class:`dict` mapping element numbers to midpoint + coordinates, with each value in the map having shape + ``(ambient_dim, nmidpoints)``. The ordering of the midpoints + follows their ordering in the tesselation. + """ + raise NotImplementedError(type(meg).__name__) + + +@singledispatch +def get_group_tesselated_nodes(meg: MeshElementGroup, tesselation, elements): + """Compute the nodes of the child elements according to the tesselation. + + :arg group: An instance of :class:`meshmode.mesh.MeshElementGroup`. + :arg tesselation: a :class:`TesselationInfo`. + :arg elements: A list of (group-relative) element numbers. + + :return: A :class:`dict` mapping element numbers to node + coordinates, with each value in the map having shape + ``(ambient_dim, nchildren, nunit_nodes)``. + The ordering of the child nodes follows the ordering + of ``tesselation.children.`` + """ + raise NotImplementedError(type(meg).__name__) + + +@singledispatch +def get_group_tesselation_info(meg: MeshElementGroup): + """ + :returns: a :class:`TesselationInfo` for the element group *meg*. + """ + raise NotImplementedError(type(meg).__name__) + +# }}} + + +# {{{ helpers + +def _midpoint_tuples(a, b): def midpoint(x, y): d, r = divmod(x + y, 2) if r: @@ -104,9 +151,7 @@ def midpoint(x, y): return tuple(midpoint(ai, bi) for ai, bi in zip(a, b)) -# {{{ resampling - -def get_ref_midpoints(shape, ref_vertices): +def _get_ref_midpoints(shape, ref_vertices): r"""The reference element is considered to be, e.g. for a 2 simplex:: F @@ -131,111 +176,83 @@ def get_ref_midpoints(shape, ref_vertices): ] return [rv for rv in ref_vertices if rv not in orig_vertices] +# }}} -def get_group_midpoints(group, tesselation, elements): - """Compute the midpoints of the vertices of the specified elements. - :arg group: an instance of :class:`meshmode.mesh.MeshElementGroup`. - :arg tesselation: a :class:`TesselationInfo`. - :arg elements: a list of (group-relative) element numbers. +# {{{ modepy.shape tesselation and resampling - :return: A :class:`dict` mapping element numbers to midpoint - coordinates, with each value in the map having shape - ``(ambient_dim, nmidpoints)``. The ordering of the midpoints - follows their ordering in the tesselation. - """ - from meshmode.mesh import _ModepyElementGroup - if not isinstance(group, _ModepyElementGroup): - raise TypeError(f"groups of type '{type(group.mesh_el_group).__name__}'" - " are not supported.") - - shape = group._modepy_shape - space = mp.space_for_shape(shape, group.order) +@get_group_midpoints.register(_ModepyElementGroup) +def _(meg: _ModepyElementGroup, tesselation, elements): + shape = meg._modepy_shape + space = meg._modepy_space # get midpoints in reference coordinates - midpoints = -1 + np.array(get_ref_midpoints(shape, tesselation.ref_vertices)) + midpoints = -1 + np.array(_get_ref_midpoints(shape, tesselation.ref_vertices)) # resample midpoints to ambient coordinates resampling_mat = mp.resampling_matrix( mp.basis_for_space(space, shape).functions, midpoints.T, - group.unit_nodes) + meg.unit_nodes) resampled_midpoints = np.einsum("mu,deu->edm", - resampling_mat, group.nodes[:, elements]) + resampling_mat, meg.nodes[:, elements]) return dict(zip(elements, resampled_midpoints)) -def get_group_tesselated_nodes(group, tesselation, elements): - """Compute the nodes of the child elements according to the tesselation. - - :arg group: An instance of :class:`meshmode.mesh.MeshElementGroup`. - :arg tesselation: a :class:`TesselationInfo`. - :arg elements: A list of (group-relative) element numbers. - - :return: A :class:`dict` mapping element numbers to node - coordinates, with each value in the map having shape - ``(ambient_dim, nchildren, nunit_nodes)``. - The ordering of the child nodes follows the ordering - of ``tesselation.children.`` - """ - from meshmode.mesh import _ModepyElementGroup - if not isinstance(group, _ModepyElementGroup): - raise TypeError(f"groups of type '{type(group.mesh_el_group).__name__}'" - " are not supported.") - - shape = group._modepy_shape - space = mp.space_for_shape(shape, group.order) +@get_group_tesselated_nodes.register(_ModepyElementGroup) +def _(meg: _ModepyElementGroup, tesselation, elements): + shape = meg._modepy_shape + space = meg._modepy_space # get child unit node coordinates. from meshmode.mesh.refinement.utils import map_unit_nodes_to_children child_unit_nodes = np.hstack(list( - map_unit_nodes_to_children(group, group.unit_nodes, tesselation) + map_unit_nodes_to_children(meg, meg.unit_nodes, tesselation) )) # resample child nodes to ambient coordinates resampling_mat = mp.resampling_matrix( mp.basis_for_space(space, shape).functions, child_unit_nodes, - group.unit_nodes) + meg.unit_nodes) resampled_unit_nodes = np.einsum("cu,deu->edc", - resampling_mat, group.nodes[:, elements]) + resampling_mat, meg.nodes[:, elements]) - ambient_dim = len(group.nodes) - nunit_nodes = len(group.unit_nodes[0]) + ambient_dim = len(meg.nodes) + nunit_nodes = len(meg.unit_nodes[0]) return { el: resampled_unit_nodes[iel].reshape((ambient_dim, -1, nunit_nodes)) for iel, el in enumerate(elements) } -# }}} - -# {{{ tesselation +@get_group_tesselation_info.register(_ModepyElementGroup) +def _(meg: _ModepyElementGroup): + shape = meg._modepy_shape + space = type(meg._modepy_space)(meg.dim, 2) -def get_shape_tesselation_info(shape): - space = mp.space_for_shape(shape, 2) ref_vertices = mp.node_tuples_for_space(space) ref_vertices_to_index = {rv: i for i, rv in enumerate(ref_vertices)} from pytools import add_tuples - space = mp.space_for_shape(shape, 1) + space = type(meg._modepy_space)(meg.dim, 1) orig_vertices = tuple([ add_tuples(vt, vt) for vt in mp.node_tuples_for_space(space) ]) orig_vertex_indices = [ref_vertices_to_index[vt] for vt in orig_vertices] - midpoints = get_ref_midpoints(shape, ref_vertices) + midpoints = _get_ref_midpoints(shape, ref_vertices) midpoint_indices = [ref_vertices_to_index[mp] for mp in midpoints] midpoint_to_vertex_pairs = { midpoint: (i, j) for i, ivt in enumerate(orig_vertices) for j, jvt in enumerate(orig_vertices) - for midpoint in [midpoint_tuples(ivt, jvt)] + for midpoint in [_midpoint_tuples(ivt, jvt)] if i < j and midpoint in midpoints } # ensure order matches the one in midpoint_indices @@ -249,17 +266,6 @@ def get_shape_tesselation_info(shape): midpoint_vertex_pairs=midpoint_vertex_pairs, ) - -@memoize -def get_group_tesselation_info(group_type, dim): - from meshmode.mesh import _ModepyElementGroup - if issubclass(group_type, _ModepyElementGroup): - shape = group_type._modepy_shape_cls(dim) - return get_shape_tesselation_info(shape) - else: - raise NotImplementedError( - "bisection for element groups of type {group_type.__name__}") - # }}} # vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index 1e31f1ab8..bfcca2ad7 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -48,7 +48,7 @@ def map_unit_nodes_to_children(meg: MeshElementGroup, raise NotImplementedError(type(meg).__name__) -@singledispatch.register(SimplexElementGroup) +@map_unit_nodes_to_children.register(SimplexElementGroup) def _(meg: SimplexElementGroup, unit_nodes, tess_info): ref_vertices = np.array(tess_info.ref_vertices, dtype=np.float).T assert len(unit_nodes.shape) == 2 @@ -61,23 +61,18 @@ def _(meg: SimplexElementGroup, unit_nodes, tess_info): yield basis.dot((unit_nodes + 1.0) / 2.0) + origin - 1.0 -@singledispatch.register(TensorProductElementGroup) +@map_unit_nodes_to_children.register(TensorProductElementGroup) def _(meg: TensorProductElementGroup, unit_nodes, tess_info): ref_vertices = np.array(tess_info.ref_vertices, dtype=np.float).T assert len(unit_nodes.shape) == 2 - import modepy as mp - from dataclass import replace - space = replace(meg._modepy_space, order=1) - basis_indices = [ - i for i, node in mp.node_tuples_for_space(space) - if sum(node) == 1 - ] + # NOTE: nodes indices in the unit hypercube that form the `e_i` basis + basis_indices = 2**np.arange(meg.dim) for child in tess_info.children: child_arr = np.array(child) - origin = ref_vertices[:, child[0]].reshape(-1, 1) - basis = ref_vertices[:, child[basis_indices]] - origin + origin = ref_vertices[:, child_arr[0]].reshape(-1, 1) + basis = ref_vertices[:, child_arr[basis_indices]] - origin # mapped nodes are on [0, 2], so we subtract 1 to get it to [-1, 1] yield basis.dot((unit_nodes + 1.0) / 2.0) + origin - 1.0 From 8cbdc82e65b65ba6f8507b09f0bf5cfacec84702 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 16 Dec 2020 15:45:29 -0600 Subject: [PATCH 11/20] rename tesselation to tess_info --- .../discretization/connection/refinement.py | 6 ++-- meshmode/mesh/refinement/__init__.py | 18 +++++------ meshmode/mesh/refinement/no_adjacency.py | 20 ++++++------ meshmode/mesh/refinement/tesselate.py | 32 +++++++++---------- 4 files changed, 38 insertions(+), 38 deletions(-) diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index df31c8f78..4c802af20 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -69,8 +69,8 @@ def _build_interpolation_batches_for_group( """ from meshmode.discretization.connection.direct import InterpolationBatch - num_children = len(record.tesselation.children) \ - if record.tesselation else 0 + num_children = len(record.tess_info.children) \ + if record.tess_info else 0 from_bins = [[] for i in range(1 + num_children)] to_bins = [[] for i in range(1 + num_children)] for elt_idx, refinement_result in enumerate(record.element_mapping): @@ -91,7 +91,7 @@ def _build_interpolation_batches_for_group( from meshmode.mesh.refinement.utils import map_unit_nodes_to_children mapped_unit_nodes = map_unit_nodes_to_children( - meg, fine_unit_nodes, record.tesselation) + meg, fine_unit_nodes, record.tess_info) from itertools import chain for from_bin, to_bin, unit_nodes in zip( diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index a9ed8050a..2f363757a 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -55,14 +55,14 @@ def get_vertex_pair_to_midpoint_order(dim): )) @staticmethod - def get_midpoints(group, tesselation, elements): + def get_midpoints(group, tess_info, elements): import meshmode.mesh.refinement.tesselate as tess - return tess.get_group_midpoints(group, tesselation, elements) + return tess.get_group_midpoints(group, tess_info, elements) @staticmethod - def get_tesselated_nodes(group, tesselation, elements): + def get_tesselated_nodes(group, tess_info, elements): import meshmode.mesh.refinement.tesselate as tess - return tess.get_group_tesselated_nodes(group, tesselation, elements) + return tess.get_group_tesselated_nodes(group, tess_info, elements) def tesselate_simplex(dim): @@ -631,7 +631,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr iel_base = grp.element_nr_base # List of lists mapping element number to new element number(s). element_mapping = [] - tesselation = None + tess_info = None # {{{ get midpoint coordinates for vertices @@ -644,7 +644,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr midpoints_to_find.append(iel_grp) if not resampler: resampler = SimplexResampler() - tesselation = TesselationInfo( + tess_info = TesselationInfo( children=self.simplex_result[grp.dim], ref_vertices=self.simplex_node_tuples[grp.dim]) else: @@ -653,7 +653,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr if midpoints_to_find: midpoints = resampler.get_midpoints( - grp, tesselation, midpoints_to_find) + grp, tess_info, midpoints_to_find) midpoint_order = resampler.get_vertex_pair_to_midpoint_order(grp.dim) del midpoints_to_find @@ -732,7 +732,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr self.group_refinement_records.append( GroupRefinementRecord( - tesselation=tesselation, + tess_info=tess_info, element_mapping=element_mapping) ) @@ -793,7 +793,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr if is_simplex: resampler = SimplexResampler() new_nodes = resampler.get_tesselated_nodes( - prev_group, refinement_record.tesselation, to_resample) + prev_group, refinement_record.tess_info, to_resample) else: raise NotImplementedError( "unimplemented: node resampling for non simplex elements") diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index dda338b08..6d89e8746 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -89,7 +89,7 @@ def refine(self, refine_flags): get_group_tesselated_nodes) for igrp, group in enumerate(mesh.groups): - tesselation = get_group_tesselation_info(group) + tess_info = get_group_tesselation_info(group) # {{{ compute counts and index arrays @@ -97,7 +97,7 @@ def refine(self, refine_flags): group.element_nr_base: group.element_nr_base+group.nelements] - nchildren = len(tesselation.children) + nchildren = len(tess_info.children) nchild_elements = np.ones(group.nelements, dtype=mesh.element_id_dtype) nchild_elements[grp_flags] = nchildren @@ -116,7 +116,7 @@ def refine(self, refine_flags): from meshmode.mesh.refinement.tesselate import GroupRefinementRecord group_refinement_records.append( GroupRefinementRecord( - tesselation=tesselation, + tess_info=tess_info, element_mapping=[ list(range( child_el_indices[iel], @@ -127,7 +127,7 @@ def refine(self, refine_flags): if perform_vertex_updates: midpoints = get_group_midpoints( - group, tesselation, refining_el_old_indices) + group, tess_info, refining_el_old_indices) new_vertex_indices = np.empty( (new_nelements, group.vertex_indices.shape[1]), @@ -141,17 +141,17 @@ def refine(self, refine_flags): for old_iel in refining_el_old_indices: new_iel_base = child_el_indices[old_iel] - refining_vertices = np.empty(len(tesselation.ref_vertices), + refining_vertices = np.empty(len(tess_info.ref_vertices), dtype=mesh.vertex_id_dtype) refining_vertices.fill(-17) # carry over old vertices - refining_vertices[tesselation.orig_vertex_indices] = \ + refining_vertices[tess_info.orig_vertex_indices] = \ group.vertex_indices[old_iel] for imidpoint, (iref_midpoint, (v1, v2)) in enumerate(zip( - tesselation.midpoint_indices, - tesselation.midpoint_vertex_pairs)): + tess_info.midpoint_indices, + tess_info.midpoint_vertex_pairs)): global_v1 = group.vertex_indices[old_iel, v1] global_v2 = group.vertex_indices[old_iel, v2] @@ -175,7 +175,7 @@ def refine(self, refine_flags): assert (refining_vertices >= 0).all() new_vertex_indices[new_iel_base:new_iel_base+nchildren] = \ - refining_vertices[tesselation.children] + refining_vertices[tess_info.children] assert (new_vertex_indices >= 0).all() else: @@ -195,7 +195,7 @@ def refine(self, refine_flags): new_nodes[:, unrefined_el_new_indices] = group.nodes[:, ~grp_flags] tesselated_nodes = get_group_tesselated_nodes( - group, tesselation, refining_el_old_indices) + group, tess_info, refining_el_old_indices) for old_iel in refining_el_old_indices: new_iel_base = child_el_indices[old_iel] diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index 43d85c1df..6ef7f6478 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -45,7 +45,7 @@ class TesselationInfo: """ .. attribute:: children - A tesselation of the reference element, given here by + A tess_info of the reference element, given here by :attr:`ref_vertices`. .. attribute:: ref_vertices @@ -81,9 +81,9 @@ class TesselationInfo: @dataclass(frozen=True) class GroupRefinementRecord: """ - .. attribute:: tesselation + .. attribute:: tess_info - A :class:`TesselationInfo` that describes the tesselation of the + A :class:`TesselationInfo` that describes the tess_info of the element group. .. attribute:: element_mapping @@ -91,39 +91,39 @@ class GroupRefinementRecord: A mapping from the original elements to the refined child elements. """ - tesselation: TesselationInfo + tess_info: TesselationInfo element_mapping: List[List[int]] @singledispatch -def get_group_midpoints(meg: MeshElementGroup, tesselation, elements): +def get_group_midpoints(meg: MeshElementGroup, tess_info, elements): """Compute the midpoints of the vertices of the specified elements. :arg group: an instance of :class:`meshmode.mesh.MeshElementGroup`. - :arg tesselation: a :class:`TesselationInfo`. + :arg tess_info: a :class:`TesselationInfo`. :arg elements: a list of (group-relative) element numbers. :return: A :class:`dict` mapping element numbers to midpoint coordinates, with each value in the map having shape ``(ambient_dim, nmidpoints)``. The ordering of the midpoints - follows their ordering in the tesselation. + follows their ordering in the tess_info. """ raise NotImplementedError(type(meg).__name__) @singledispatch -def get_group_tesselated_nodes(meg: MeshElementGroup, tesselation, elements): - """Compute the nodes of the child elements according to the tesselation. +def get_group_tesselated_nodes(meg: MeshElementGroup, tess_info, elements): + """Compute the nodes of the child elements according to the tess_info. :arg group: An instance of :class:`meshmode.mesh.MeshElementGroup`. - :arg tesselation: a :class:`TesselationInfo`. + :arg tess_info: a :class:`TesselationInfo`. :arg elements: A list of (group-relative) element numbers. :return: A :class:`dict` mapping element numbers to node coordinates, with each value in the map having shape ``(ambient_dim, nchildren, nunit_nodes)``. The ordering of the child nodes follows the ordering - of ``tesselation.children.`` + of ``tess_info.children.`` """ raise NotImplementedError(type(meg).__name__) @@ -179,15 +179,15 @@ def _get_ref_midpoints(shape, ref_vertices): # }}} -# {{{ modepy.shape tesselation and resampling +# {{{ modepy.shape tess_info and resampling @get_group_midpoints.register(_ModepyElementGroup) -def _(meg: _ModepyElementGroup, tesselation, elements): +def _(meg: _ModepyElementGroup, tess_info, elements): shape = meg._modepy_shape space = meg._modepy_space # get midpoints in reference coordinates - midpoints = -1 + np.array(_get_ref_midpoints(shape, tesselation.ref_vertices)) + midpoints = -1 + np.array(_get_ref_midpoints(shape, tess_info.ref_vertices)) # resample midpoints to ambient coordinates resampling_mat = mp.resampling_matrix( @@ -202,14 +202,14 @@ def _(meg: _ModepyElementGroup, tesselation, elements): @get_group_tesselated_nodes.register(_ModepyElementGroup) -def _(meg: _ModepyElementGroup, tesselation, elements): +def _(meg: _ModepyElementGroup, tess_info, elements): shape = meg._modepy_shape space = meg._modepy_space # get child unit node coordinates. from meshmode.mesh.refinement.utils import map_unit_nodes_to_children child_unit_nodes = np.hstack(list( - map_unit_nodes_to_children(meg, meg.unit_nodes, tesselation) + map_unit_nodes_to_children(meg, meg.unit_nodes, tess_info) )) # resample child nodes to ambient coordinates From cfb3f7054dee4446f063afc2baed76ecf9918d59 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 16 Dec 2020 15:53:48 -0600 Subject: [PATCH 12/20] fix docs after rename --- meshmode/mesh/refinement/tesselate.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index 6ef7f6478..860b663a6 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -83,7 +83,7 @@ class GroupRefinementRecord: """ .. attribute:: tess_info - A :class:`TesselationInfo` that describes the tess_info of the + A :class:`TesselationInfo` that describes the tesselation of the element group. .. attribute:: element_mapping @@ -106,14 +106,14 @@ def get_group_midpoints(meg: MeshElementGroup, tess_info, elements): :return: A :class:`dict` mapping element numbers to midpoint coordinates, with each value in the map having shape ``(ambient_dim, nmidpoints)``. The ordering of the midpoints - follows their ordering in the tess_info. + follows their ordering in the tesselation. """ raise NotImplementedError(type(meg).__name__) @singledispatch def get_group_tesselated_nodes(meg: MeshElementGroup, tess_info, elements): - """Compute the nodes of the child elements according to the tess_info. + """Compute the nodes of the child elements according to the tesselation. :arg group: An instance of :class:`meshmode.mesh.MeshElementGroup`. :arg tess_info: a :class:`TesselationInfo`. @@ -179,7 +179,7 @@ def _get_ref_midpoints(shape, ref_vertices): # }}} -# {{{ modepy.shape tess_info and resampling +# {{{ modepy.shape tesselation and resampling @get_group_midpoints.register(_ModepyElementGroup) def _(meg: _ModepyElementGroup, tess_info, elements): From 23694c40780cb7e0daf64a25ca7e36cd9fca54be Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 3 Jan 2021 20:23:08 -0600 Subject: [PATCH 13/20] rename meg to fine_meg MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- meshmode/discretization/connection/refinement.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 4c802af20..828390d60 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -87,11 +87,11 @@ def _build_interpolation_batches_for_group( to_bin.append(child_idx) fine_unit_nodes = fine_discr_group.unit_nodes - meg = fine_discr_group.mesh_el_group + fine_meg = fine_discr_group.mesh_el_group from meshmode.mesh.refinement.utils import map_unit_nodes_to_children mapped_unit_nodes = map_unit_nodes_to_children( - meg, fine_unit_nodes, record.tess_info) + fine_meg, fine_unit_nodes, record.tess_info) from itertools import chain for from_bin, to_bin, unit_nodes in zip( From e46d4f2bba7e05abb93c0677f7cc4e8978c9dfe0 Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 3 Jan 2021 20:23:32 -0600 Subject: [PATCH 14/20] fix type in docs MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- meshmode/mesh/refinement/tesselate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index 860b663a6..25445d826 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -152,7 +152,7 @@ def midpoint(x, y): def _get_ref_midpoints(shape, ref_vertices): - r"""The reference element is considered to be, e.g. for a 2 simplex:: + r"""The reference element is considered to be, e.g. for a 2-simplex:: F | \ From ae730008acbffc88c88b2f76c681fa4a66d106ad Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 3 Jan 2021 20:24:34 -0600 Subject: [PATCH 15/20] add fixme for element_mapping MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- meshmode/mesh/refinement/tesselate.py | 1 + 1 file changed, 1 insertion(+) diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index 25445d826..c6e62c721 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -92,6 +92,7 @@ class GroupRefinementRecord: """ tess_info: TesselationInfo + # FIXME: This should really be a CSR data structure. element_mapping: List[List[int]] From ff165844c2c34cd86c23089f11d15b390794e11d Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 3 Jan 2021 20:24:57 -0600 Subject: [PATCH 16/20] improve docs for TesselationInfo MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- meshmode/mesh/refinement/tesselate.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index c6e62c721..d2bad78e6 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -41,8 +41,9 @@ # {{{ interface @dataclass(frozen=True) -class TesselationInfo: - """ +class ElementTesselationInfo: + """Describes how one element is split into multiple child elements. + .. attribute:: children A tess_info of the reference element, given here by From 06f0974faea3d5a13c19eb7e25a9115e5c041fae Mon Sep 17 00:00:00 2001 From: Alex Fikl Date: Sun, 3 Jan 2021 20:25:25 -0600 Subject: [PATCH 17/20] improve docs for TesselationInfo MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Andreas Klöckner --- meshmode/mesh/refinement/tesselate.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index d2bad78e6..efd3a382b 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -84,8 +84,8 @@ class GroupRefinementRecord: """ .. attribute:: tess_info - A :class:`TesselationInfo` that describes the tesselation of the - element group. + An instance of :class:`ElementTesselationInfo` that describes the + tesselation of a single element into multiple child elements. .. attribute:: element_mapping From da7596257c4f8afe338d5eab96490e6680b484a0 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 Jan 2021 20:28:48 -0600 Subject: [PATCH 18/20] improve docs for ElementTesselationInfo --- meshmode/mesh/refinement/tesselate.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index efd3a382b..d4b2c2551 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -43,11 +43,11 @@ @dataclass(frozen=True) class ElementTesselationInfo: """Describes how one element is split into multiple child elements. - + .. attribute:: children - A tess_info of the reference element, given here by - :attr:`ref_vertices`. + An array of shape ``(nchildren, nvertices)`` containing the vertices + of each child element the reference element was split into. .. attribute:: ref_vertices From ce36d6a040a14ed40dbb61685b2f0be2bd9cf89d Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 Jan 2021 20:58:31 -0600 Subject: [PATCH 19/20] rename to el_tess_info --- meshmode/mesh/refinement/__init__.py | 20 +++++++++--------- meshmode/mesh/refinement/no_adjacency.py | 20 +++++++++--------- meshmode/mesh/refinement/tesselate.py | 26 ++++++++++++------------ meshmode/mesh/refinement/utils.py | 17 ++++++++-------- 4 files changed, 42 insertions(+), 41 deletions(-) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index 2f363757a..c3d258d20 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -27,7 +27,7 @@ from meshmode.mesh.refinement.no_adjacency import RefinerWithoutAdjacency from meshmode.mesh.refinement.tesselate import \ - TesselationInfo, GroupRefinementRecord + ElementTesselationInfo, GroupRefinementRecord import logging logger = logging.getLogger(__name__) @@ -55,14 +55,14 @@ def get_vertex_pair_to_midpoint_order(dim): )) @staticmethod - def get_midpoints(group, tess_info, elements): + def get_midpoints(group, el_tess_info, elements): import meshmode.mesh.refinement.tesselate as tess - return tess.get_group_midpoints(group, tess_info, elements) + return tess.get_group_midpoints(group, el_tess_info, elements) @staticmethod - def get_tesselated_nodes(group, tess_info, elements): + def get_tesselated_nodes(group, el_tess_info, elements): import meshmode.mesh.refinement.tesselate as tess - return tess.get_group_tesselated_nodes(group, tess_info, elements) + return tess.get_group_tesselated_nodes(group, el_tess_info, elements) def tesselate_simplex(dim): @@ -631,7 +631,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr iel_base = grp.element_nr_base # List of lists mapping element number to new element number(s). element_mapping = [] - tess_info = None + el_tess_info = None # {{{ get midpoint coordinates for vertices @@ -644,7 +644,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr midpoints_to_find.append(iel_grp) if not resampler: resampler = SimplexResampler() - tess_info = TesselationInfo( + el_tess_info = ElementTesselationInfo( children=self.simplex_result[grp.dim], ref_vertices=self.simplex_node_tuples[grp.dim]) else: @@ -653,7 +653,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr if midpoints_to_find: midpoints = resampler.get_midpoints( - grp, tess_info, midpoints_to_find) + grp, el_tess_info, midpoints_to_find) midpoint_order = resampler.get_vertex_pair_to_midpoint_order(grp.dim) del midpoints_to_find @@ -732,7 +732,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr self.group_refinement_records.append( GroupRefinementRecord( - tess_info=tess_info, + el_tess_info=el_tess_info, element_mapping=element_mapping) ) @@ -793,7 +793,7 @@ def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_gr if is_simplex: resampler = SimplexResampler() new_nodes = resampler.get_tesselated_nodes( - prev_group, refinement_record.tess_info, to_resample) + prev_group, refinement_record.el_tess_info, to_resample) else: raise NotImplementedError( "unimplemented: node resampling for non simplex elements") diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index 6d89e8746..cc9caa37e 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -89,7 +89,7 @@ def refine(self, refine_flags): get_group_tesselated_nodes) for igrp, group in enumerate(mesh.groups): - tess_info = get_group_tesselation_info(group) + el_tess_info = get_group_tesselation_info(group) # {{{ compute counts and index arrays @@ -97,7 +97,7 @@ def refine(self, refine_flags): group.element_nr_base: group.element_nr_base+group.nelements] - nchildren = len(tess_info.children) + nchildren = len(el_tess_info.children) nchild_elements = np.ones(group.nelements, dtype=mesh.element_id_dtype) nchild_elements[grp_flags] = nchildren @@ -116,7 +116,7 @@ def refine(self, refine_flags): from meshmode.mesh.refinement.tesselate import GroupRefinementRecord group_refinement_records.append( GroupRefinementRecord( - tess_info=tess_info, + el_tess_info=el_tess_info, element_mapping=[ list(range( child_el_indices[iel], @@ -127,7 +127,7 @@ def refine(self, refine_flags): if perform_vertex_updates: midpoints = get_group_midpoints( - group, tess_info, refining_el_old_indices) + group, el_tess_info, refining_el_old_indices) new_vertex_indices = np.empty( (new_nelements, group.vertex_indices.shape[1]), @@ -141,17 +141,17 @@ def refine(self, refine_flags): for old_iel in refining_el_old_indices: new_iel_base = child_el_indices[old_iel] - refining_vertices = np.empty(len(tess_info.ref_vertices), + refining_vertices = np.empty(len(el_tess_info.ref_vertices), dtype=mesh.vertex_id_dtype) refining_vertices.fill(-17) # carry over old vertices - refining_vertices[tess_info.orig_vertex_indices] = \ + refining_vertices[el_tess_info.orig_vertex_indices] = \ group.vertex_indices[old_iel] for imidpoint, (iref_midpoint, (v1, v2)) in enumerate(zip( - tess_info.midpoint_indices, - tess_info.midpoint_vertex_pairs)): + el_tess_info.midpoint_indices, + el_tess_info.midpoint_vertex_pairs)): global_v1 = group.vertex_indices[old_iel, v1] global_v2 = group.vertex_indices[old_iel, v2] @@ -175,7 +175,7 @@ def refine(self, refine_flags): assert (refining_vertices >= 0).all() new_vertex_indices[new_iel_base:new_iel_base+nchildren] = \ - refining_vertices[tess_info.children] + refining_vertices[el_tess_info.children] assert (new_vertex_indices >= 0).all() else: @@ -195,7 +195,7 @@ def refine(self, refine_flags): new_nodes[:, unrefined_el_new_indices] = group.nodes[:, ~grp_flags] tesselated_nodes = get_group_tesselated_nodes( - group, tess_info, refining_el_old_indices) + group, el_tess_info, refining_el_old_indices) for old_iel in refining_el_old_indices: new_iel_base = child_el_indices[old_iel] diff --git a/meshmode/mesh/refinement/tesselate.py b/meshmode/mesh/refinement/tesselate.py index d4b2c2551..13b50ee33 100644 --- a/meshmode/mesh/refinement/tesselate.py +++ b/meshmode/mesh/refinement/tesselate.py @@ -82,7 +82,7 @@ class ElementTesselationInfo: @dataclass(frozen=True) class GroupRefinementRecord: """ - .. attribute:: tess_info + .. attribute:: el_tess_info An instance of :class:`ElementTesselationInfo` that describes the tesselation of a single element into multiple child elements. @@ -92,17 +92,17 @@ class GroupRefinementRecord: A mapping from the original elements to the refined child elements. """ - tess_info: TesselationInfo + el_tess_info: ElementTesselationInfo # FIXME: This should really be a CSR data structure. element_mapping: List[List[int]] @singledispatch -def get_group_midpoints(meg: MeshElementGroup, tess_info, elements): +def get_group_midpoints(meg: MeshElementGroup, el_tess_info, elements): """Compute the midpoints of the vertices of the specified elements. :arg group: an instance of :class:`meshmode.mesh.MeshElementGroup`. - :arg tess_info: a :class:`TesselationInfo`. + :arg el_tess_info: a :class:`ElementTesselationInfo`. :arg elements: a list of (group-relative) element numbers. :return: A :class:`dict` mapping element numbers to midpoint @@ -114,18 +114,18 @@ def get_group_midpoints(meg: MeshElementGroup, tess_info, elements): @singledispatch -def get_group_tesselated_nodes(meg: MeshElementGroup, tess_info, elements): +def get_group_tesselated_nodes(meg: MeshElementGroup, el_tess_info, elements): """Compute the nodes of the child elements according to the tesselation. :arg group: An instance of :class:`meshmode.mesh.MeshElementGroup`. - :arg tess_info: a :class:`TesselationInfo`. + :arg el_tess_info: a :class:`ElementTesselationInfo`. :arg elements: A list of (group-relative) element numbers. :return: A :class:`dict` mapping element numbers to node coordinates, with each value in the map having shape ``(ambient_dim, nchildren, nunit_nodes)``. The ordering of the child nodes follows the ordering - of ``tess_info.children.`` + of ``el_tess_info.children.`` """ raise NotImplementedError(type(meg).__name__) @@ -133,7 +133,7 @@ def get_group_tesselated_nodes(meg: MeshElementGroup, tess_info, elements): @singledispatch def get_group_tesselation_info(meg: MeshElementGroup): """ - :returns: a :class:`TesselationInfo` for the element group *meg*. + :returns: a :class:`ElementTesselationInfo` for the element group *meg*. """ raise NotImplementedError(type(meg).__name__) @@ -184,12 +184,12 @@ def _get_ref_midpoints(shape, ref_vertices): # {{{ modepy.shape tesselation and resampling @get_group_midpoints.register(_ModepyElementGroup) -def _(meg: _ModepyElementGroup, tess_info, elements): +def _(meg: _ModepyElementGroup, el_tess_info, elements): shape = meg._modepy_shape space = meg._modepy_space # get midpoints in reference coordinates - midpoints = -1 + np.array(_get_ref_midpoints(shape, tess_info.ref_vertices)) + midpoints = -1 + np.array(_get_ref_midpoints(shape, el_tess_info.ref_vertices)) # resample midpoints to ambient coordinates resampling_mat = mp.resampling_matrix( @@ -204,14 +204,14 @@ def _(meg: _ModepyElementGroup, tess_info, elements): @get_group_tesselated_nodes.register(_ModepyElementGroup) -def _(meg: _ModepyElementGroup, tess_info, elements): +def _(meg: _ModepyElementGroup, el_tess_info, elements): shape = meg._modepy_shape space = meg._modepy_space # get child unit node coordinates. from meshmode.mesh.refinement.utils import map_unit_nodes_to_children child_unit_nodes = np.hstack(list( - map_unit_nodes_to_children(meg, meg.unit_nodes, tess_info) + map_unit_nodes_to_children(meg, meg.unit_nodes, el_tess_info) )) # resample child nodes to ambient coordinates @@ -260,7 +260,7 @@ def _(meg: _ModepyElementGroup): # ensure order matches the one in midpoint_indices midpoint_vertex_pairs = [midpoint_to_vertex_pairs[m] for m in midpoints] - return TesselationInfo( + return ElementTesselationInfo( ref_vertices=ref_vertices, children=np.array(mp.submesh_for_shape(shape, ref_vertices)), orig_vertex_indices=np.array(orig_vertex_indices), diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index bfcca2ad7..4ccde721f 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -37,11 +37,12 @@ @singledispatch def map_unit_nodes_to_children(meg: MeshElementGroup, - unit_nodes, tess_info) -> np.ndarray: + unit_nodes, el_tess_info) -> np.ndarray: """ :arg unit_nodes: an :class:`~numpy.ndarray` of unit nodes on the element type described by *meg*. - :arg tess_info: a :class:`~meshmode.mesh.refinement.tesselate.TesselationInfo`. + :arg el_tess_info: a + :class:`~meshmode.mesh.refinement.tesselate.ElementTesselationInfo`. :returns: an :class:`~numpy.ndarray` of mapped unit nodes for each child in the tesselation. """ @@ -49,11 +50,11 @@ def map_unit_nodes_to_children(meg: MeshElementGroup, @map_unit_nodes_to_children.register(SimplexElementGroup) -def _(meg: SimplexElementGroup, unit_nodes, tess_info): - ref_vertices = np.array(tess_info.ref_vertices, dtype=np.float).T +def _(meg: SimplexElementGroup, unit_nodes, el_tess_info): + ref_vertices = np.array(el_tess_info.ref_vertices, dtype=np.float).T assert len(unit_nodes.shape) == 2 - for child in tess_info.children: + for child in el_tess_info.children: origin = ref_vertices[:, child[0]].reshape(-1, 1) basis = ref_vertices[:, child[1:]] - origin @@ -62,14 +63,14 @@ def _(meg: SimplexElementGroup, unit_nodes, tess_info): @map_unit_nodes_to_children.register(TensorProductElementGroup) -def _(meg: TensorProductElementGroup, unit_nodes, tess_info): - ref_vertices = np.array(tess_info.ref_vertices, dtype=np.float).T +def _(meg: TensorProductElementGroup, unit_nodes, el_tess_info): + ref_vertices = np.array(el_tess_info.ref_vertices, dtype=np.float).T assert len(unit_nodes.shape) == 2 # NOTE: nodes indices in the unit hypercube that form the `e_i` basis basis_indices = 2**np.arange(meg.dim) - for child in tess_info.children: + for child in el_tess_info.children: child_arr = np.array(child) origin = ref_vertices[:, child_arr[0]].reshape(-1, 1) basis = ref_vertices[:, child_arr[basis_indices]] - origin From 8d8914d8e8d462faa5c05214d18b7a7bdb0dd3c5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 3 Jan 2021 21:24:21 -0600 Subject: [PATCH 20/20] rename some missed tess_info --- meshmode/discretization/connection/refinement.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 828390d60..b4c6f6c1e 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -69,8 +69,8 @@ def _build_interpolation_batches_for_group( """ from meshmode.discretization.connection.direct import InterpolationBatch - num_children = len(record.tess_info.children) \ - if record.tess_info else 0 + num_children = len(record.el_tess_info.children) \ + if record.el_tess_info else 0 from_bins = [[] for i in range(1 + num_children)] to_bins = [[] for i in range(1 + num_children)] for elt_idx, refinement_result in enumerate(record.element_mapping): @@ -91,7 +91,7 @@ def _build_interpolation_batches_for_group( from meshmode.mesh.refinement.utils import map_unit_nodes_to_children mapped_unit_nodes = map_unit_nodes_to_children( - fine_meg, fine_unit_nodes, record.tess_info) + fine_meg, fine_unit_nodes, record.el_tess_info) from itertools import chain for from_bin, to_bin, unit_nodes in zip(