diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 5ef91fe7e41dc464061e73438796afa174756560..ff009fce6efe7833e7a8bbd43c6d33745e722704 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -31,39 +31,6 @@ 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 - -# }}} - - # {{{ Build interpolation batches for group def _build_interpolation_batches_for_group( @@ -123,7 +90,8 @@ def _build_interpolation_batches_for_group( to_bin.append(child_idx) fine_unit_nodes = fine_discr_group.unit_nodes - mapped_unit_nodes = _map_unit_nodes_to_children( + from meshmode.mesh.refinement.utils import map_unit_nodes_to_children + mapped_unit_nodes = map_unit_nodes_to_children( fine_unit_nodes, record.tesselation) from itertools import chain diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index ce2a5f236d0ec524a411c78c87b7904abd423db6..ff5d8581d5d74291e48ecafe5936548b290e0960 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -547,7 +547,7 @@ def generate_warped_rect_mesh(dim, order, n): 0.05*np.cos(10*x[0]) + 1.3*x[1] + np.sin(x[1])) if len(x) == 3: - result[2] = x[2] + np.sin(x[0]) + result[2] = x[2] + np.sin(x[0] / 2) / 2 return result from meshmode.mesh.processing import map_mesh diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index 16639442bdbaf10d6dd2def5c87bd8ff3254c25a..3c39dde0a25ba950ae94cf063a3d87a56c01d9b9 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -26,8 +26,6 @@ import itertools from six.moves import range from pytools import RecordWithoutPickling -from meshmode.mesh.generation import make_group_from_vertices - class TreeRayNode(object): """Describes a ray as a tree, this class represents each node in this tree @@ -217,31 +215,6 @@ class Refiner(object): return self.previous_mesh def get_current_mesh(self): - - from meshmode.mesh import Mesh - #return Mesh(vertices, [grp], nodal_adjacency=self.generate_nodal_adjacency(len(self.last_mesh.groups[group].vertex_indices) \ - # + count*3)) - groups = [] - grpn = 0 - for grp in self.last_mesh.groups: - groups.append(np.empty([len(grp.vertex_indices), - len(self.last_mesh.groups[grpn].vertex_indices[0])], dtype=np.int32)) - for iel_grp in range(grp.nelements): - for i in range(0, len(grp.vertex_indices[iel_grp])): - groups[grpn][iel_grp][i] = grp.vertex_indices[iel_grp][i] - grpn += 1 - grp = [] - - for grpn in range(0, len(groups)): - grp.append(make_group_from_vertices(self.last_mesh.vertices, groups[grpn], 4)) - self.last_mesh = Mesh( - self.last_mesh.vertices, grp, - nodal_adjacency=self.generate_nodal_adjacency( - len(self.last_mesh.groups[0].vertex_indices), - len(self.last_mesh.vertices[0])), - vertex_id_dtype=self.last_mesh.vertex_id_dtype, - element_id_dtype=self.last_mesh.element_id_dtype) - return self.last_mesh def get_leaves(self, cur_node): @@ -590,15 +563,42 @@ class Refiner(object): element_mapping = [] tesselation = None + # {{{ get midpoint coordinates for vertices + + midpoints_to_find = [] + resampler = None + for iel_grp in range(grp.nelements): + if refine_flags[iel_base + iel_grp]: + # if simplex + 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 = self._Tesselation( + self.simplex_result[grp.dim], + self.simplex_node_tuples[grp.dim]) + else: + raise NotImplementedError("unimplemented: midpoint finding" + "for non simplex elements") + + if midpoints_to_find: + midpoints = resampler.get_midpoints( + grp, tesselation, midpoints_to_find) + midpoint_order = resampler.get_vertex_pair_to_midpoint_order(grp.dim) + + del midpoints_to_find + + # }}} + 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 simplex if len(vertex_indices) == grp.dim + 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]) @@ -617,18 +617,15 @@ class Refiner(object): vertex_pair2 = (max_index, vertices_index) self.pair_map[vertex_pair1] = cur_node.left self.pair_map[vertex_pair2] = cur_node.right - for k in range(len(self.last_mesh.vertices)): - vertices[k, vertices_index] = \ - (self.last_mesh.vertices[k, vertex_indices[i]] + - self.last_mesh.vertices[k, vertex_indices[j]]) / 2.0 + midpoint_idx = midpoint_order[(i, j)] + vertices[:, vertices_index] = \ + midpoints[iel_grp][:, midpoint_idx] midpoint_vertices.append(vertices_index) vertices_index += 1 else: cur_midpoint = cur_node.midpoint midpoint_vertices.append(cur_midpoint) - # }}} - #generate new rays cur_dim = len(grp.vertex_indices[0])-1 for i in range(len(midpoint_vertices)): @@ -653,12 +650,12 @@ class Refiner(object): for j in range(len(self.simplex_result[cur_dim][i])): 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: + else: #quadrilateral + raise NotImplementedError("unimplemented: " + "support for quad elements") # node_tuple_to_coord = {} # for node_index, node_tuple in enumerate(self.index_to_node_tuple[cur_dim]): # node_tuple_to_coord[node_tuple] = grp.vertex_indices[iel_grp][node_index] @@ -699,15 +696,64 @@ class Refiner(object): #check_adjacent_elements(groups, new_hanging_vertex_element, nelements_in_grp) self.hanging_vertex_element = new_hanging_vertex_element - grp = [] - for grpn in range(0, len(groups)): - grp.append(make_group_from_vertices(vertices, groups[grpn], 4)) + + # {{{ make new groups + + new_mesh_el_groups = [] + + for refinement_record, group, prev_group in zip( + self.group_refinement_records, groups, self.last_mesh.groups): + is_simplex = len(prev_group.vertex_indices[0]) == prev_group.dim + 1 + ambient_dim = len(prev_group.nodes) + nelements = len(group) + nunit_nodes = len(prev_group.unit_nodes[0]) + + nodes = np.empty( + (ambient_dim, nelements, nunit_nodes), + dtype=prev_group.nodes.dtype) + + element_mapping = refinement_record.element_mapping + + to_resample = [elem for elem in range(len(element_mapping)) + if len(element_mapping[elem]) > 1] + + 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) + else: + raise NotImplementedError( + "unimplemented: node resampling for non simplex elements") + + for elem, mapped_elems in enumerate(element_mapping): + if len(mapped_elems) == 1: + # No resampling required, just copy over + nodes[:, mapped_elems[0]] = prev_group.nodes[:, elem] + n = nodes[:, mapped_elems[0]] + else: + nodes[:, mapped_elems] = new_nodes[elem] + + if is_simplex: + new_mesh_el_groups.append( + type(prev_group)( + order=prev_group.order, + vertex_indices=group, + nodes=nodes, + unit_nodes=prev_group.unit_nodes)) + else: + raise NotImplementedError("unimplemented: support for creating" + "non simplex element groups") + + # }}} from meshmode.mesh import Mesh self.previous_mesh = self.last_mesh self.last_mesh = Mesh( - vertices, grp, + vertices, new_mesh_el_groups, nodal_adjacency=self.generate_nodal_adjacency( totalnelements, nvertices, groups), vertex_id_dtype=self.last_mesh.vertex_id_dtype, diff --git a/meshmode/mesh/refinement/resampler.py b/meshmode/mesh/refinement/resampler.py new file mode 100644 index 0000000000000000000000000000000000000000..cea22b3aa3f7c36b23f823e06646bb27e2723110 --- /dev/null +++ b/meshmode/mesh/refinement/resampler.py @@ -0,0 +1,134 @@ +from __future__ import division, absolute_import, print_function + +__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 modepy as mp + +import logging +logger = logging.getLogger(__name__) + + +# {{{ resampling simplex points for refinement + +# NOTE: Class internal to refiner: do not make documentation public. +class SimplexResampler(object): + """ + 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()`. + """ + + def get_vertex_pair_to_midpoint_order(self, 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) + )) + + def get_midpoints(self, 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`) + """ + 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), + midpoints.T, + group.unit_nodes) + + resamp_midpoints = np.einsum("mu,deu->edm", + resamp_mat, + group.nodes[:, elements]) + + return dict(zip(elements, resamp_midpoints)) + + def get_tesselated_nodes(self, 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.`` + """ + 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), + child_unit_nodes, + group.unit_nodes) + + resamp_unit_nodes = np.einsum("cu,deu->edc", + resamp_mat, + group.nodes[:, elements]) + + ambient_dim = len(group.nodes) + nunit_nodes = len(group.unit_nodes[0]) + + return dict((elem, + resamp_unit_nodes[ielem].reshape( + (ambient_dim, -1, nunit_nodes))) + for ielem, elem in enumerate(elements)) + +# }}} + + +# vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index 751142aa6d67ff7717b017ebe83f13a18e216ef3..08b71f8755f4f270756c904d8d5a592e753cb30c 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -29,8 +29,42 @@ 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: diff --git a/test/test_refinement.py b/test/test_refinement.py index bed61eb94c900db0f78b526cfe43ffababe4d825..8098766d446874b08e676d304d5ac73b88e1163a 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -47,10 +47,10 @@ logger = logging.getLogger(__name__) from functools import partial -def gen_blob_mesh(h=0.2): +def gen_blob_mesh(h=0.2, order=1): from meshmode.mesh.io import generate_gmsh, FileSource return generate_gmsh( - FileSource("blob-2d.step"), 2, order=1, + FileSource("blob-2d.step"), 2, order=order, force_ambient_dim=2, other_options=[ "-string", "Mesh.CharacteristicLengthMax = %s;" % h] @@ -152,16 +152,21 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): ("warp", 2, [4, 5, 6]), ("warp", 3, [4, 5, 6]), ]) +@pytest.mark.parametrize("mesh_order", [1, 5]) @pytest.mark.parametrize("refine_flags", [ # FIXME: slow - # uniform_refine_flags, + #uniform_refine_flags, partial(random_refine_flags, 0.4) ]) def test_refinement_connection( - ctx_getter, group_factory, mesh_name, dim, mesh_pars, refine_flags): + ctx_getter, group_factory, mesh_name, dim, mesh_pars, mesh_order, + refine_flags, plot_mesh=False): from random import seed seed(13) + # Discretization order + order = 5 + cl_ctx = ctx_getter() queue = cl.CommandQueue(cl_ctx) @@ -172,8 +177,6 @@ def test_refinement_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) @@ -185,14 +188,17 @@ def test_refinement_connection( assert dim == 1 h = 1 / mesh_par mesh = make_curve_mesh( - partial(ellipse, 1), np.linspace(0, 1, mesh_par + 1), order=1) + partial(ellipse, 1), np.linspace(0, 1, mesh_par + 1), + order=mesh_order) elif mesh_name == "blob": + if mesh_order == 5: + pytest.xfail("https://gitlab.tiker.net/inducer/meshmode/issues/2") assert dim == 2 h = mesh_par - mesh = gen_blob_mesh(h) + mesh = gen_blob_mesh(h, mesh_order) elif mesh_name == "warp": from meshmode.mesh.generation import generate_warped_rect_mesh - mesh = generate_warped_rect_mesh(dim, order=1, n=mesh_par) + mesh = generate_warped_rect_mesh(dim, order=mesh_order, n=mesh_par) h = 1/mesh_par else: raise ValueError("mesh_name not recognized") @@ -202,7 +208,9 @@ def test_refinement_connection( discr = Discretization(cl_ctx, mesh, group_factory(order)) refiner = Refiner(mesh) - refiner.refine(refine_flags(mesh)) + flags = refine_flags(mesh) + refiner.refine(flags) + connection = make_refinement_connection( refiner, discr, group_factory(order)) check_connection(connection) @@ -215,6 +223,18 @@ def test_refinement_connection( f_interp = connection(queue, f_coarse).with_queue(queue) f_true = f(x_fine).with_queue(queue) + if plot_mesh: + import matplotlib.pyplot as plt + x = x.get(queue) + err = np.array(np.log10( + 1e-16 + np.abs((f_interp - f_true).get(queue))), dtype=float) + import matplotlib.cm as cm + cmap = cm.ScalarMappable(cmap=cm.jet) + cmap.set_array(err) + plt.scatter(x[0], x[1], c=cmap.to_rgba(err), s=20, cmap=cmap) + plt.colorbar(cmap) + plt.show() + import numpy.linalg as la err = la.norm((f_interp - f_true).get(queue), np.inf) eoc_rec.add_data_point(h, err)