diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 8ea924efd1d7a5235747745d8e99f54a5a51a612..ba15835834a1dd7afbc53bc8c35ff258d7778e42 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -36,6 +36,9 @@ from meshmode.discretization.connection.face import ( make_face_restriction, make_face_to_all_faces_embedding) from meshmode.discretization.connection.opposite_face import \ make_opposite_face_connection +from meshmode.discretization.connection.refinement import \ + make_refinement_connection + import logging logger = logging.getLogger(__name__) @@ -47,7 +50,8 @@ __all__ = [ "FRESTR_INTERIOR_FACES", "FRESTR_ALL_FACES", "make_face_restriction", "make_face_to_all_faces_embedding", - "make_opposite_face_connection" + "make_opposite_face_connection", + "make_refinement_connection" ] __doc__ = """ @@ -62,6 +66,8 @@ __doc__ = """ .. autofunction:: make_opposite_face_connection +.. autofunction:: make_refinement_connection + Implementation details ^^^^^^^^^^^^^^^^^^^^^^ @@ -434,11 +440,4 @@ def check_connection(connection): # }}} -# {{{ refinement connection - -def make_refinement_connection(refiner, coarse_discr): - pass - -# }}} - # vim: foldmethod=marker diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py new file mode 100644 index 0000000000000000000000000000000000000000..75604ec05dd945cde45e0e6cdd55724f88083c3f --- /dev/null +++ b/meshmode/discretization/connection/refinement.py @@ -0,0 +1,195 @@ +# -*- coding: utf-8 -*- +from __future__ import division, print_function, absolute_import + +__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. +""" + +import numpy as np +import pyopencl as cl +import pyopencl.array # noqa + +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) + + 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 + +# }}} + + +# {{{ Build interpolation batches for group + +def _build_interpolation_batches_for_group( + queue, group_idx, coarse_discr_group, fine_discr_group, record): + """ + To map between discretizations, we sort each of the fine mesh + elements into an interpolation batch. Which batch they go + into is determined by where the refined unit nodes live + relative to the coarse reference element. + + For instance, consider the following refinement: + + ______ ______ + |\ | |\ e| + | \ | |d\ | + | \ | |__\ | + | \ | => |\c|\ | + | \ | |a\|b\ | + | \| | | \| + ‾‾‾‾‾‾ ‾‾‾‾‾‾ + + Here, the discretization unit nodes for elements a,b,c,d,e + will each have different positions relative to the reference + element, so each element gets its own batch. On the other + hand, for + + ______ ______ + |\ | |\ f|\e| + | \ | |d\ |g\| + | \ | |__\|__| + | \ | => |\c|\ | + | \ | |a\|b\h| + | \| | | \| + ‾‾‾‾‾‾ ‾‾‾‾‾‾ + + the pairs {a,e}, {b,f}, {c,g}, {d,h} can share interpolation + batches because their unit nodes are mapped from the same part + of the reference element. + """ + from meshmode.discretization.connection import InterpolationBatch + + num_children = len(record.tesselation.children) \ + if record.tesselation 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): + if len(refinement_result) == 1: + # Not refined -> interpolates to self + from_bins[0].append(elt_idx) + to_bins[0].append(refinement_result[0]) + else: + assert len(refinement_result) == num_children + # Refined -> interpolates to children + for from_bin, to_bin, child_idx in zip( + from_bins[1:], to_bins[1:], refinement_result): + from_bin.append(elt_idx) + to_bin.append(child_idx) + + fine_unit_nodes = fine_discr_group.unit_nodes + mapped_unit_nodes = _map_unit_nodes_to_children( + fine_unit_nodes, record.tesselation) + + from itertools import chain + for from_bin, to_bin, unit_nodes in zip( + from_bins, to_bins, + chain([fine_unit_nodes], mapped_unit_nodes)): + if not from_bin: + continue + yield InterpolationBatch( + from_group_index=group_idx, + from_element_indices=cl.array.to_device(queue, np.asarray(from_bin)), + to_element_indices=cl.array.to_device(queue, np.asarray(to_bin)), + result_unit_nodes=unit_nodes, + to_element_face=None) + +# }}} + + +def make_refinement_connection(refiner, coarse_discr, group_factory): + """Return a + :class:`meshmode.discretization.connection.DiscretizationConnection` + connecting `coarse_discr` to a discretization on the fine mesh. + + :arg refiner: An instance of + :class:`meshmode.mesh.refinement.Refiner` + + :arg coarse_discr: An instance of + :class:`meshmode.mesh.discretization.Discretization` associated + with the mesh given to the refiner + + :arg group_factory: An instance of + :class:`meshmode.mesh.discretization.ElementGroupFactory`. Used for + discretizing the fine mesh. + """ + from meshmode.discretization.connection import ( + DiscretizationConnectionElementGroup, DiscretizationConnection) + + coarse_mesh = refiner.get_previous_mesh() + fine_mesh = refiner.last_mesh + assert coarse_discr.mesh is coarse_mesh + + from meshmode.discretization import Discretization + fine_discr = Discretization( + coarse_discr.cl_context, + fine_mesh, + group_factory, + real_dtype=coarse_discr.real_dtype) + + logger.info("building refinement connection: start") + + groups = [] + with cl.CommandQueue(fine_discr.cl_context) as queue: + for group_idx, (coarse_discr_group, fine_discr_group, record) in \ + enumerate(zip(coarse_discr.groups, fine_discr.groups, + refiner.group_refinement_records)): + groups.append( + DiscretizationConnectionElementGroup( + list(_build_interpolation_batches_for_group( + queue, group_idx, coarse_discr_group, + fine_discr_group, record)))) + + logger.info("building refinement connection: done") + + return DiscretizationConnection( + from_discr=coarse_discr, + to_discr=fine_discr, + groups=groups, + is_surjective=True) + + +# vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index 0430ada9929bb270e230bf183d324141e890804c..2154d45195412c6c012a54ace9b8fbcb07561055 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -24,6 +24,7 @@ THE SOFTWARE. import numpy as np import itertools from six.moves import range +from pytools import RecordWithoutPickling from meshmode.mesh.generation import make_group_from_vertices @@ -57,18 +58,31 @@ class TreeRayNode(object): class Refiner(object): + class _Tesselation(RecordWithoutPickling): + + def __init__(self, children, ref_vertices): + RecordWithoutPickling.__init__(self, + ref_vertices=ref_vertices, children=children) + + class _GroupRefinementRecord(RecordWithoutPickling): + + def __init__(self, tesselation, element_mapping): + RecordWithoutPickling.__init__(self, + tesselation=tesselation, element_mapping=element_mapping) + # {{{ constructor def __init__(self, mesh): from meshmode.mesh.tesselate import tesselatetet, tesselatetri self.lazy = False self.seen_tuple = {} + self.group_refinement_records = [] tri_node_tuples, tri_result = tesselatetri() tet_node_tuples, tet_result = tesselatetet() - print(tri_node_tuples, tri_result) #quadrilateral_node_tuples = [ #print tri_result, tet_result self.simplex_node_tuples = [None, None, tri_node_tuples, tet_node_tuples] + # Dimension-parameterized tesselations for refinement self.simplex_result = [None, None, tri_result, tet_result] #print tri_node_tuples, tri_result #self.simplex_node_tuples, self.simplex_result = tesselatetet() @@ -89,6 +103,8 @@ class Refiner(object): for i in range(nvertices): self.hanging_vertex_element.append([]) + # Fill pair_map. + # Add adjacency information to each TreeRayNode. for grp in mesh.groups: iel_base = grp.element_nr_base for iel_grp in range(grp.nelements): @@ -195,6 +211,9 @@ class Refiner(object): self.last_mesh.nelements - self.get_refine_base_index(), np.bool) + def get_previous_mesh(self): + return self.previous_mesh + def get_current_mesh(self): from meshmode.mesh import Mesh @@ -417,6 +436,8 @@ class Refiner(object): # for next_vertices in next_vertices_list: # remove_element_from_connectivity(next_vertices, new_hanging_vertex_elements, to_remove) + # {{{ Add element to connectivity + def add_element_to_connectivity(vertices, new_hanging_vertex_elements, to_add): if len(vertices) == 2: min_vertex = min(vertices[0], vertices[1]) @@ -502,14 +523,22 @@ class Refiner(object): # return # add_element_to_connectivity(next_element_rays, new_hanging_vertex_elements, to_add) + # }}} + + # {{{ Add hanging vertex element + def add_hanging_vertex_el(v_index, el): assert not (v_index == 37 and el == 48) new_hanging_vertex_element[v_index].append(el) + # }}} + # def remove_ray_el(ray, el): # ray.remove(el) + # {{{ Check adjacent elements + def check_adjacent_elements(groups, new_hanging_vertex_elements, nelements_in_grp): for grp in groups: iel_base = 0 @@ -529,6 +558,8 @@ class Refiner(object): assert((iel_base+iel_grp) in new_hanging_vertex_elements[cur_node.left_vertex]) assert((iel_base+iel_grp) in new_hanging_vertex_elements[cur_node.right_vertex]) + # }}} + for i in range(len(self.last_mesh.vertices)): for j in range(len(self.last_mesh.vertices[i])): vertices[i,j] = self.last_mesh.vertices[i,j] @@ -545,14 +576,24 @@ class Refiner(object): grpn = 0 vertices_index = len(self.last_mesh.vertices[0]) nelements_in_grp = grp.nelements - for grp in self.last_mesh.groups: + del self.group_refinement_records[:] + + for grp_idx, grp in enumerate(self.last_mesh.groups): iel_base = grp.element_nr_base + # List of lists mapping element number to new element number(s). + element_mapping = [] + tesselation = None + for iel_grp in range(grp.nelements): + element_mapping.append([iel_grp]) if refine_flags[iel_base+iel_grp]: midpoint_vertices = [] vertex_indices = grp.vertex_indices[iel_grp] #if simplex if len(grp.vertex_indices[iel_grp]) == len(self.last_mesh.vertices)+1: + + # {{{ Get midpoints for all pairs of vertices + for i in range(len(vertex_indices)): for j in range(i+1, len(vertex_indices)): min_index = min(vertex_indices[i], vertex_indices[j]) @@ -581,6 +622,7 @@ class Refiner(object): cur_midpoint = cur_node.midpoint midpoint_vertices.append(cur_midpoint) + # }}} #generate new rays cur_dim = len(grp.vertex_indices[0])-1 @@ -598,13 +640,16 @@ class Refiner(object): for midpoint_index, midpoint_tuple in enumerate(self.index_to_midpoint_tuple[cur_dim]): node_tuple_to_coord[midpoint_tuple] = midpoint_vertices[midpoint_index] for i in range(len(self.simplex_result[cur_dim])): + if i == 0: + iel = iel_grp + else: + iel = nelements_in_grp + i - 1 + element_mapping[-1].append(iel) for j in range(len(self.simplex_result[cur_dim][i])): - if i == 0: - groups[grpn][iel_grp][j] = \ - node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]] - else: - groups[grpn][nelements_in_grp+i-1][j] = \ - node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]] + groups[grpn][iel][j] = \ + node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]] + tesselation = self._Tesselation( + self.simplex_result[cur_dim], self.simplex_node_tuples[cur_dim]) nelements_in_grp += len(self.simplex_result[cur_dim])-1 #assuming quad otherwise #else: @@ -615,6 +660,8 @@ class Refiner(object): # def generate_all_tuples(cur_list): # if len(cur_list[len(cur_list)-1]) + self.group_refinement_records.append( + self._GroupRefinementRecord(tesselation, element_mapping)) #clear connectivity data for grp in self.last_mesh.groups: @@ -653,6 +700,7 @@ class Refiner(object): from meshmode.mesh import Mesh + self.previous_mesh = self.last_mesh self.last_mesh = Mesh( vertices, grp, nodal_adjacency=self.generate_nodal_adjacency( diff --git a/test/test_refinement.py b/test/test_refinement.py index aa1ab68a725c73c8618e32fa7c7a8dd518a4420b..9c36b9750ebe87a899a8dd5c5f17f4191ed1b0b4 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -23,6 +23,8 @@ THE SOFTWARE. """ import pytest +import pyopencl as cl +import pyopencl.clmath # noqa import numpy as np from pyopencl.tools import ( # noqa @@ -33,19 +35,25 @@ from meshmode.mesh.generation import ( # noqa from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry from meshmode.mesh.refinement import Refiner +from meshmode.discretization.poly_element import ( + InterpolatoryQuadratureSimplexGroupFactory, + PolynomialWarpAndBlendGroupFactory, + PolynomialEquidistantGroupFactory, +) + import logging logger = logging.getLogger(__name__) from functools import partial -def gen_blob_mesh(): +def gen_blob_mesh(h=0.2): from meshmode.mesh.io import generate_gmsh, FileSource return generate_gmsh( FileSource("blob-2d.step"), 2, order=1, force_ambient_dim=2, other_options=[ - "-string", "Mesh.CharacteristicLengthMax = %s;" % 0.2] + "-string", "Mesh.CharacteristicLengthMax = %s;" % h] ) @@ -124,6 +132,84 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): check_nodal_adj_against_geometry(mesh) +@pytest.mark.parametrize("group_factory", [ + InterpolatoryQuadratureSimplexGroupFactory, + PolynomialWarpAndBlendGroupFactory, + PolynomialEquidistantGroupFactory + ]) +@pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ + ("blob", 2, [1e-1, 8e-2, 5e-2]), + ("warp", 2, [4, 5, 6]), + ("warp", 3, [4, 5, 6]), +]) +@pytest.mark.parametrize("refine_flags", [ + # FIXME: slow + # uniform_refine_flags, + partial(random_refine_flags, 0.4) +]) +def test_refinement_connection( + ctx_getter, group_factory, mesh_name, dim, mesh_pars, refine_flags): + from random import seed + seed(13) + + cl_ctx = ctx_getter() + queue = cl.CommandQueue(cl_ctx) + + from meshmode.discretization import Discretization + from meshmode.discretization.connection import ( + make_refinement_connection, check_connection) + + from pytools.convergence import EOCRecorder + eoc_rec = EOCRecorder() + + order = 5 + + def f(x): + from six.moves import reduce + return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + + for mesh_par in mesh_pars: + # {{{ get mesh + + if mesh_name == "blob": + assert dim == 2 + h = mesh_par + mesh = gen_blob_mesh(h) + elif mesh_name == "warp": + from meshmode.mesh.generation import generate_warped_rect_mesh + mesh = generate_warped_rect_mesh(dim, order=1, n=mesh_par) + h = 1/mesh_par + else: + raise ValueError("mesh_name not recognized") + + # }}} + + discr = Discretization(cl_ctx, mesh, group_factory(order)) + + refiner = Refiner(mesh) + refiner.refine(refine_flags(mesh)) + connection = make_refinement_connection( + refiner, discr, group_factory(order)) + check_connection(connection) + + fine_discr = connection.to_discr + + x = discr.nodes().with_queue(queue) + x_fine = fine_discr.nodes().with_queue(queue) + f_coarse = f(x) + f_interp = connection(queue, f_coarse).with_queue(queue) + f_true = f(x_fine).with_queue(queue) + + import numpy.linalg as la + err = la.norm((f_interp - f_true).get(queue), np.inf) + eoc_rec.add_data_point(h, err) + + print(eoc_rec) + assert ( + eoc_rec.order_estimate() >= order-0.5 + or eoc_rec.max_error() < 1e-14) + + if __name__ == "__main__": import sys if len(sys.argv) > 1: