From 939673301613cb6881ae382a80706b22294cef81 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sat, 1 Jun 2019 20:50:32 -0500 Subject: [PATCH] Come up with a design for how mesh periodicity might be represented, and how to create meshes without having to compute the full connectivity --- meshmode/mesh/__init__.py | 75 ++++++++++++++++++++++++++++++++---- meshmode/mesh/periodicity.py | 67 ++++++++++++++++++++++++++++++++ meshmode/mesh/tools.py | 18 +++++++++ 3 files changed, 153 insertions(+), 7 deletions(-) create mode 100644 meshmode/mesh/periodicity.py diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 6d6191aa..be5af1b2 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -618,6 +618,17 @@ class Mesh(Record): (Note that element groups are not necessarily contiguous like the figure may suggest.) + .. attribute:: facial_boundary_adjacency_groups + + A list of either None or :class:`FacialAdjacencyGroup` instances, one + per element group. If *None*, boundary information for the + corresponding element group is not yet known. + + Once :attr:`facial_adjacency_groups` is computed, corresponds to the + value for the key ``None`` in each group's connectivity dictionary. + + .. versionadded:: 2019.1 + .. attribute:: boundary_tags A tuple of boundary tag identifiers. :class:`BTAG_ALL` and @@ -649,6 +660,8 @@ class Mesh(Record): skip_element_orientation_test=False, nodal_adjacency=None, facial_adjacency_groups=None, + facial_boundary_adjacency_groups=None, + periodicity=None, boundary_tags=None, vertex_id_dtype=np.int32, element_id_dtype=np.int32, @@ -669,19 +682,29 @@ class Mesh(Record): will be deduced from vertex adjacency. *False*, in which case this information will be marked unavailable (such as if there are hanging nodes in the geometry, so that vertex adjacency does not convey - the full picture), and references to - :attr:`element_neighbors_starts` and :attr:`element_neighbors` - will result in exceptions. Lastly, a tuple + the full picture). Lastly, a tuple :class:`NodalAdjacency` object. :arg facial_adjacency_groups: One of three options: *None*, in which case this information will be deduced from vertex adjacency. *False*, in which case this information will be marked unavailable (such as if there are hanging nodes in the geometry, so that vertex adjacency does not convey - the full picture), and references to - :attr:`element_neighbors_starts` and :attr:`element_neighbors` - will result in exceptions. Lastly, a data structure as described in + the full picture). Lastly, a data structure as described in :attr:`facial_adjacency_groups` may be passed. + :arg facial_boundary_adjacency_groups: One of three options: + *None*, in which case this information + will be deduced from vertex adjacency. *False*, in which case + this information will be marked unavailable (such as if there are + hanging nodes in the geometry, so that vertex adjacency does not convey + the full picture). Lastly, a data structure as described in + :attr:`facial_boundary_adjacency_groups` may be passed. + :arg periodicity: *None* or a list of + :class:`meshmode.mesh.periodicity.Periodicity` instances, + describing connectivity across a periodic mesh boundary. + + .. versionchanged:: 2019.1 + + *periodicity* and *facial_boundary_adjacency_groups* were added. """ el_nr = 0 @@ -725,6 +748,8 @@ class Mesh(Record): nodal_adjacency = False if facial_adjacency_groups is None: facial_adjacency_groups = False + if facial_boundary_adjacency_groups is None: + facial_boundary_adjacency_groups = False if nodal_adjacency is not False and nodal_adjacency is not None: if not isinstance(nodal_adjacency, NodalAdjacency): @@ -740,6 +765,9 @@ class Mesh(Record): self, vertices=vertices, groups=new_groups, _nodal_adjacency=nodal_adjacency, _facial_adjacency_groups=facial_adjacency_groups, + _facial_boundary_adjacency_groups=( + facial_boundary_adjacency_groups), + _periodicity=periodicity, boundary_tags=boundary_tags, btag_to_index=btag_to_index, vertex_id_dtype=np.dtype(vertex_id_dtype), @@ -805,6 +833,9 @@ class Mesh(Record): set_if_not_present("boundary_tags") set_if_not_present("nodal_adjacency", "_nodal_adjacency") set_if_not_present("facial_adjacency_groups", "_facial_adjacency_groups") + set_if_not_present("facial_boundary_adjacency_groups", + "_facial_boundary_adjacency_groups") + set_if_not_present("periodicity", "_periodicity") set_if_not_present("vertex_id_dtype") set_if_not_present("element_id_dtype") @@ -870,6 +901,30 @@ class Mesh(Record): return self._facial_adjacency_groups + @property + def facial_boundary_adjacency_groups(self): + def compute_from_facial_adj(): + assert self._facial_boundary_adjacency_groups is None + return [fagrps.get(None, None) + for fagrps in self._facial_adjacency_groups] + + if self._facial_adjacency_groups is not None: + return compute_from_facial_adj() + + if self._facial_boundary_adjacency_groups is False: + from meshmode import DataUnavailable + raise DataUnavailable("facial_boundary_adjacency_groups") + + elif self._facial_boundary_adjacency_groups is None: + if not self.is_conforming: + from meshmode import DataUnavailable + raise DataUnavailable("facial_boundary_adjacency_groups can only " + "be computed for known-conforming meshes") + + return compute_from_facial_adj() + + return self._facial_adjacency_groups + def boundary_tag_bit(self, boundary_tag): try: return 1 << self.btag_to_index[boundary_tag] @@ -887,6 +942,10 @@ class Mesh(Record): == other._nodal_adjacency) and (self._facial_adjacency_groups == other._facial_adjacency_groups) + and (self._facial_boundary_adjacency_groups + == other._facial_boundary_adjacency_groups) + and (self._periodicity + == other._periodicity) and self.boundary_tags == other.boundary_tags and self.is_conforming == other.is_conforming) @@ -1046,6 +1105,7 @@ def _compute_facial_adjacency_from_vertices(groups, boundary_tags, del igrp # {{{ build facial_adjacency_groups data structure, still empty + from meshmode.mesh import FacialAdjacencyGroup, BTAG_ALL, BTAG_REALLY_ALL facial_adjacency_groups = [] @@ -1254,7 +1314,8 @@ def as_python(mesh, function_name="make_mesh"): cg(" boundary_tags=[%s]," % btags_str) cg(" is_conforming=%s)" % repr(mesh.is_conforming)) - # FIXME: Handle facial adjacency, boundary tags + # FIXME: Handle facial adjacency, boundary tags, + # facial_boundary_adjacency_groups return cg.get() diff --git a/meshmode/mesh/periodicity.py b/meshmode/mesh/periodicity.py new file mode 100644 index 00000000..6b311d80 --- /dev/null +++ b/meshmode/mesh/periodicity.py @@ -0,0 +1,67 @@ +from __future__ import division, absolute_import + +__copyright__ = "Copyright (C) 2019 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. +""" + +# from six.moves import range +# import six +# +# import numpy as np +# import modepy as mp +# import numpy.linalg as la +from pytools import Record + + +class Periodicity(Record): + """From the perspective of the user, this should be considered an opaque + data structure, to be passed to the periodicities argument of the + :class:`meshmode.mesh.Mesh` constructor. + """ + + __doc_private__ = """Periodicity is considered to have two sides, called + "A" and "B" for definiteness. + + .. attribute:: a_to_b_affine_map + + .. attribute:: b_to_a_vertex_map + + A :class:`dict` mapping vertices on the "B" side of the periodicity + to the "A" side of the periodicity. + """ + + def __init__(self, a_to_b_affine_map, b_to_a_vertex_map): + self.a_to_b_affine_map = a_to_b_affine_map + self.b_to_a_vertex_map = b_to_a_vertex_map + + +# {{{ connect_periodicity + +def connect_periodicity(mesh, boundary_tag_a, boundary_tag_b, a_to_b_affine_map): + from meshmode.mesh import make_vertex_lookup_tree + vlookup = make_vertex_lookup_tree(mesh) # noqa + + # FIXME + raise NotImplementedError() + +# }}} + +# vim: foldmethod=marker diff --git a/meshmode/mesh/tools.py b/meshmode/mesh/tools.py index df1c353c..f07892ce 100644 --- a/meshmode/mesh/tools.py +++ b/meshmode/mesh/tools.py @@ -52,6 +52,24 @@ def make_element_lookup_tree(mesh, eps=1e-12): # }}} +# {{{ make_vertex_lookup_tree + +def make_vertex_lookup_tree(mesh, eps=1e-12): + from meshmode.mesh.processing import find_bounding_box + bbox_min, bbox_max = find_bounding_box(mesh) + bbox_min -= eps + bbox_max += eps + + tree = SpatialBinaryTreeBucket(bbox_min, bbox_max) + + for ivert, vert in enumerate(mesh.vertices.T): + tree.insert(ivert, (vert-eps, vert+eps)) + + return tree + +# }}} + + # {{{ nd_quad_submesh def nd_quad_submesh(node_tuples): -- GitLab