diff --git a/examples/refinement-playground.py b/examples/refinement-playground.py new file mode 100644 index 0000000000000000000000000000000000000000..a29c0db8f687440f3493015cb256d30926f6b4d8 --- /dev/null +++ b/examples/refinement-playground.py @@ -0,0 +1,305 @@ +from __future__ import division, print_function + +from six.moves import range +import numpy as np # noqa +import pyopencl as cl +import random +import os +import logging +order = 1 + +from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry +from meshmode.mesh.refinement import Refiner + + +#construct vertex vertex_index +def remove_if_exists(name): + from errno import ENOENT + try: + os.remove(name) + except OSError as e: + if e.errno == ENOENT: + pass + else: + raise + + +def linear_func(vert): + csum = 0 + for i in vert: + csum += i + #print csum + return csum + + +def sine_func(vert): + #print vert + import math + res = 1 + for i in vert: + res *= math.sin(i*2.0*math.pi) + #print res + return abs(res) * 0.2 + 0.2 + + +def quadratic_func(vert): + csum = 0 + for i in vert: + csum += i * i + return csum * 1.5 + + +def get_function_flags(mesh, function): + from math import sqrt + flags = np.zeros(len(mesh.groups[0].vertex_indices)) + for grp in mesh.groups: + for iel_grp in range(grp.nelements): + vertex_indices = grp.vertex_indices[iel_grp] + max_edge_len = 0 + for i in range(len(vertex_indices)): + for j in range(i+1, len(vertex_indices)): + edge_len = 0 + for k in range(len(mesh.vertices)): + edge_len += ( + (mesh.vertices[k, vertex_indices[i]] + - mesh.vertices[k, vertex_indices[j]]) + * (mesh.vertices[k, vertex_indices[i]] + - mesh.vertices[k, vertex_indices[j]])) + edge_len = sqrt(edge_len) + max_edge_len = max(max_edge_len, edge_len) + #print(edge_lens[0], mesh.vertices[0, vertex_indices[i]], mesh.vertices[1, vertex_indices[i]], mesh.vertices[2, vertex_indices[i]]) # noqa + centroid = [0] * len(mesh.vertices) + for j in range(len(mesh.vertices)): + centroid[j] += mesh.vertices[j, vertex_indices[i]] + for i in range(len(mesh.vertices)): + centroid[i] /= len(vertex_indices) + val = function(centroid) + if max_edge_len > val: + flags[iel_grp] = True + return flags + + +def get_corner_flags(mesh): + flags = np.zeros(len(mesh.groups[0].vertex_indices)) + for grp in mesh.groups: + for iel_grp in range(grp.nelements): + is_corner_el = False + vertex_indices = grp.vertex_indices[iel_grp] + for i in range(len(vertex_indices)): + cur_vertex_corner = True + for j in range(len(mesh.vertices)): + print(iel_grp, i, mesh.vertices[j, vertex_indices[i]]) + if mesh.vertices[j, vertex_indices[i]] != 0.0: + cur_vertex_corner = False + if cur_vertex_corner: + is_corner_el = True + break + if is_corner_el: + print(iel_grp) + flags[iel_grp] = True + return flags + + +def get_random_flags(mesh): + flags = np.zeros(len(mesh.groups[0].vertex_indices)) + for i in range(0, len(flags)): + flags[i] = random.randint(0, 1) + return flags + + +def refine_and_generate_chart_function(mesh, filename, function): + from time import clock + cl_ctx = cl.create_some_context() + queue = cl.CommandQueue(cl_ctx) + print("NELEMENTS: ", mesh.nelements) + #print mesh + for i in range(len(mesh.groups[0].vertex_indices[0])): + for k in range(len(mesh.vertices)): + print(mesh.vertices[k, i]) + + #check_nodal_adj_against_geometry(mesh); + r = Refiner(mesh) + #random.seed(0) + #times = 3 + num_elements = [] + time_t = [] + #nelements = mesh.nelements + while True: + print("NELS:", mesh.nelements) + #flags = get_corner_flags(mesh) + flags = get_function_flags(mesh, function) + nels = 0 + for i in flags: + if i: + nels += 1 + if nels == 0: + break + print("LKJASLFKJALKASF:", nels) + num_elements.append(nels) + #flags = get_corner_flags(mesh) + beg = clock() + mesh = r.refine(flags) + end = clock() + time_taken = end - beg + time_t.append(time_taken) + #if nelements == mesh.nelements: + #break + #nelements = mesh.nelements + #from meshmode.mesh.visualization import draw_2d_mesh + #draw_2d_mesh(mesh, True, True, True, fill=None) + #import matplotlib.pyplot as pt + #pt.show() + + #poss_flags = np.zeros(len(mesh.groups[0].vertex_indices)) + #for i in range(0, len(flags)): + # poss_flags[i] = flags[i] + #for i in range(len(flags), len(poss_flags)): + # poss_flags[i] = 1 + + import matplotlib.pyplot as pt + pt.xlabel('Number of elements being refined') + pt.ylabel('Time taken') + pt.plot(num_elements, time_t, "o") + pt.savefig(filename, format='pdf') + pt.clf() + print('DONE REFINING') + ''' + flags = np.zeros(len(mesh.groups[0].vertex_indices)) + flags[0] = 1 + flags[1] = 1 + mesh = r.refine(flags) + flags = np.zeros(len(mesh.groups[0].vertex_indices)) + flags[0] = 1 + flags[1] = 1 + flags[2] = 1 + mesh = r.refine(flags) + ''' + #check_nodal_adj_against_geometry(mesh) + #r.print_rays(70) + #r.print_rays(117) + #r.print_hanging_elements(10) + #r.print_hanging_elements(117) + #r.print_hanging_elements(757) + #from meshmode.mesh.visualization import draw_2d_mesh + #draw_2d_mesh(mesh, False, False, False, fill=None) + #import matplotlib.pyplot as pt + #pt.show() + + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + PolynomialWarpAndBlendGroupFactory + discr = Discretization( + cl_ctx, mesh, PolynomialWarpAndBlendGroupFactory(order)) + from meshmode.discretization.visualization import make_visualizer + vis = make_visualizer(queue, discr, order) + remove_if_exists("connectivity2.vtu") + remove_if_exists("geometry2.vtu") + vis.write_vtk_file("geometry2.vtu", [ + ("f", discr.nodes()[0]), + ]) + + from meshmode.discretization.visualization import \ + write_nodal_adjacency_vtk_file + + write_nodal_adjacency_vtk_file("connectivity2.vtu", + mesh) + + +def main2(): + from meshmode.mesh.generation import ( # noqa + generate_icosphere, generate_icosahedron, + generate_torus, generate_regular_rect_mesh, + generate_box_mesh) +# mesh = generate_icosphere(1, order=order) + #mesh = generate_icosahedron(1, order=order) + #mesh = generate_torus(3, 1, order=order) + #mesh = generate_regular_rect_mesh() + #mesh = generate_box_mesh(3*(np.linspace(0, 3, 5),)) + #mesh = generate_box_mesh(3*(np.linspace(0, 1, 3),)) + mesh = generate_box_mesh(3*(np.linspace(0, 1, 5),)) + refine_and_generate_chart_function(mesh, "plot.pdf", sine_func) + + +def all_refine(num_mesh, depth, fname): + from meshmode.mesh.generation import ( # noqa + generate_icosphere, generate_icosahedron, + generate_torus, generate_regular_rect_mesh, + generate_box_mesh) + import timeit + nelements = [] + runtimes = [] + for el_fact in range(2, num_mesh+2): + mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),)) + r = Refiner(mesh) + for time in range(depth): + flags = np.ones(len(mesh.groups[0].vertex_indices)) + if time < depth-1: + mesh = r.refine(flags) + else: + start = timeit.default_timer() + mesh = r.refine(flags) + stop = timeit.default_timer() + nelements.append(mesh.nelements) + runtimes.append(stop-start) + check_nodal_adj_against_geometry(mesh) + import matplotlib.pyplot as pt + pt.plot(nelements, runtimes, "o") + pt.savefig(fname) + pt.clf() + #pt.show() + + +def uniform_refine(num_mesh, fract, depth, fname): + from meshmode.mesh.generation import ( # noqa + generate_icosphere, generate_icosahedron, + generate_torus, generate_regular_rect_mesh, + generate_box_mesh) + import timeit + nelements = [] + runtimes = [] + for el_fact in range(2, num_mesh+2): + mesh = generate_box_mesh(3*(np.linspace(0, 1, el_fact),)) + r = Refiner(mesh) + all_els = list(range(mesh.nelements)) + for time in range(depth): + print("EL_FACT", el_fact, "TIME", time) + flags = np.zeros(mesh.nelements) + from random import shuffle, seed + seed(1) + shuffle(all_els) + nels_this_round = 0 + for i in range(len(all_els)): + if i / len(flags) > fract: + break + flags[all_els[i]] = 1 + nels_this_round += 1 + + if time < depth-1: + mesh = r.refine(flags) + else: + start = timeit.default_timer() + mesh = r.refine(flags) + stop = timeit.default_timer() + nelements.append(mesh.nelements) + runtimes.append(stop-start) + all_els = [] + for i in range(len(flags)): + if flags[i]: + all_els.append(i) + for i in range(len(flags), mesh.nelements): + all_els.append(i) + check_nodal_adj_against_geometry(mesh) + + import matplotlib.pyplot as pt + pt.plot(nelements, runtimes, "o") + pt.savefig(fname) + pt.clf() + + +if __name__ == "__main__": + logging.basicConfig(level="DEBUG") + print("HEREERERE") + #all_refine(3, 2, 'all_a.pdf') + all_refine(3, 3, 'all_b.pdf') + #uniform_refine(3, 0.2, 3, 'uniform_a.pdf') + #main2() diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 269ed5d14c61e3b33ef3d7e6d4354b105e9b5f50..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 ^^^^^^^^^^^^^^^^^^^^^^ @@ -103,7 +109,7 @@ class InterpolationBatch(object): .. attribute:: result_unit_nodes A :class:`numpy.ndarray` of shape - ``(from_group.dim,to_group.nelements,to_group.nunit_nodes)`` + ``(from_group.dim,to_group.nunit_nodes)`` storing the coordinates of the nodes (in unit coordinates of the *from* reference element) from which the node locations of this element should be interpolated. @@ -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/discretization/visualization.py b/meshmode/discretization/visualization.py index d68ec0c58078e68417721c8205fc3f24021d0bad..956e301533b3222db7d50916ed381222b26e416a 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -291,7 +291,7 @@ def write_nodal_adjacency_vtk_file(file_name, mesh, compressor=None,): nconnections = len(adj.neighbors) connections = np.empty((nconnections, 2), dtype=np.int32) - nb_starts = adj.neighbors_starts + nb_starts = adj.neighbor_starts for iel in range(mesh.nelements): connections[nb_starts[iel]:nb_starts[iel+1], 0] = iel diff --git a/meshmode/mesh/refinement.py b/meshmode/mesh/refinement.py deleted file mode 100644 index efe92c7197d344dd814f7966da1022201c967a12..0000000000000000000000000000000000000000 --- a/meshmode/mesh/refinement.py +++ /dev/null @@ -1,185 +0,0 @@ -from __future__ import division - -__copyright__ = "Copyright (C) 2014 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 - - -# {{{ simple, no-connectivity refinement - -def _refine_simplex_el_group(grp, refine_flags, vertex_nr_base): - from pytools import ( - generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam) - from modepy.tools import submesh - - node_tuples = gnitstam(2, grp.dim) - template = submesh(node_tuples) - - nrefined_elements = np.sum(refine_flags != 0) - - old_vertex_tuples = ( - [(0,) * grp.dim] - + [(0,)*i + (2,) + (0,)*(grp.dim-1-i) for i in range(grp.dim)]) - new_vertex_tuples = [ - vertex_tuple - for vertex_tuple in node_tuples - if vertex_tuple not in old_vertex_tuples] - - new_unit_nodes = np.array(new_vertex_tuples, dtype=grp.nodes.dtype) - 1 - - nnew_vertices = nrefined_elements * len(new_vertex_tuples) - ambient_dim = grp.nodes.shape[0] - new_vertices = np.empty((nnew_vertices, ambient_dim)) - - nnew_elements = (len(template) - 1) * nrefined_elements - new_vertex_indices = np.empty( - (nnew_elements, grp.dim+1), - dtype=grp.vertex_indices.dtype) - - new_nodes = np.empty( - (ambient_dim, nnew_elements, grp.nunit_nodes), - dtype=grp.nodes.dtype) - - inew_vertex = vertex_nr_base - inew_el = 0 - for iel in xrange(grp.nelements): - if not refine_flags[iel]: - new_nodes[:, inew_el, :] = grp.nodes[:, iel, :] - - # old vertices always keep their numbers - new_vertex_indices[inew_el, :] = grp.vertex_indices[iel, :] - continue - - el_vertex_indices = grp.vertex_indices[iel, :] - - - - - - - - - - - - - from meshmode.mesh import SimplexElementGroup - new_grp = SimplexElementGroup(grp.order, new_vertex_indices, - new_nodes, unit_nodes=grp.unit_nodes) - - return new_vertices, new_grp - - -def refine_without_connectivity(mesh, refine_flags): - vertex_chunks = [mesh.vertices] - nvertices = mesh.vertices.shape[1] - - groups = [] - - from meshmode.mesh import SimplexElementGroup, Mesh - for grp in mesh.groups: - if isinstance(grp, SimplexElementGroup): - enb = grp.element_nr_base - added_vertices, new_grp = _refine_simplex_el_group( - grp, refine_flags[enb:enb+grp.nelements], nvertices) - - vertex_chunks.append(added_vertices) - groups.append(new_grp) - else: - raise NotImplementedError("refining element group of type %s" - % type(grp).__name__) - - vertices = np.hstack(vertex_chunks) - - if all(grp.dim == 1 for grp in groups): - # For meshes that are exclusively 1D, deducing connectivity from vertex - # indices is still OK. For everyone else, not really. - - connectivity = None - else: - connectivity = False - - return Mesh(vertices, groups, element_connectivity=connectivity) - -# }}} - - -class _SplitFaceRecord(object): - """ - .. attribute:: neighboring_elements - .. attribute:: new_vertex_nr - - integer or None - """ - - -class Refiner(object): - def __init__(self, mesh): - self.last_mesh = mesh - - # Indices in last_mesh that were split in the last round of - # refinement. Only these elements may be split this time - # around. - self.last_split_elements = None - - def get_refine_base_index(self): - if self.last_split_elements is None: - return 0 - else: - return self.last_mesh.nelements - len(self.last_split_elements) - - def get_empty_refine_flags(self): - return np.zeros( - self.last_mesh.nelements - self.get_refine_base_index(), - np.bool) - - def refine(self, refine_flags): - """ - :return: a refined mesh - """ - - # multi-level mapping: - # { - # dimension of intersecting entity (0=vertex,1=edge,2=face,...) - # : - # { frozenset of end vertices : _SplitFaceRecord } - # } - split_faces = {} - - ibase = self.get_refine_base_index() - affected_group_indices = set() - - for grp in self.last_mesh.groups: - iel_base - - - - - - - - - - - -# vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py new file mode 100644 index 0000000000000000000000000000000000000000..2154d45195412c6c012a54ace9b8fbcb07561055 --- /dev/null +++ b/meshmode/mesh/refinement/__init__.py @@ -0,0 +1,853 @@ +from __future__ import division, print_function + +__copyright__ = "Copyright (C) 2014-6 Shivam Gupta" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + + +import numpy as np +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 + .. attribute:: left + Left child. + *None* if ray segment hasn't been split yet. + .. attribute:: right + Right child. + *None* if ray segment hasn't been split yet. + .. attribute:: midpoint + Vertex index of midpoint of this ray segment. + *None* if no midpoint has been assigned yet. + .. attribute:: adjacent_elements + List containing elements indices of elements adjacent + to this ray segment. + """ + def __init__(self, left_vertex, right_vertex, adjacent_elements=[]): + import copy + self.left = None + self.right = None + self.parent = None + self.left_vertex = left_vertex + self.right_vertex = right_vertex + self.midpoint = None + self.adjacent_elements = copy.deepcopy(adjacent_elements) + self.adjacent_add_diff = [] + + +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() + #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() + self.last_mesh = mesh + + # {{{ initialization + + # mapping: (vertex_1, vertex_2) -> TreeRayNode + # where vertex_i represents a vertex number + # + # Assumption: vertex_1 < vertex_2 + self.pair_map = {} + + nvertices = len(mesh.vertices[0]) + + #array containing element whose edge lies on corresponding vertex + self.hanging_vertex_element = [] + 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): + vert_indices = grp.vertex_indices[iel_grp] + for i in range(len(vert_indices)): + for j in range(i+1, len(vert_indices)): + + #minimum and maximum of the two indices for storing + #data in vertex_pair + min_index = min(vert_indices[i], vert_indices[j]) + max_index = max(vert_indices[i], vert_indices[j]) + + vertex_pair = (min_index, max_index) + #print(vertex_pair) + if vertex_pair not in self.pair_map: + self.pair_map[vertex_pair] = TreeRayNode(min_index, max_index) + self.pair_map[vertex_pair].adjacent_elements.append(iel_base+iel_grp) + elif (iel_base+iel_grp) not in self.pair_map[vertex_pair].adjacent_elements: + (self.pair_map[vertex_pair]. + adjacent_elements.append(iel_base+iel_grp)) + # }}} + + #print(vert_indices) + #generate reference tuples + self.index_to_node_tuple = [] + self.index_to_midpoint_tuple = [] + for d in range(len(vert_indices)): + dim = d + 1 + cur_index_to_node_tuple = [()] * dim + for i in range(0, dim-1): + cur_index_to_node_tuple[0] = cur_index_to_node_tuple[0] + (0,) + for i in range(1, dim): + for j in range(1, dim): + if i == j: + cur_index_to_node_tuple[i] = cur_index_to_node_tuple[i] + (2,) + else: + cur_index_to_node_tuple[i] = cur_index_to_node_tuple[i] + (0,) + cur_index_to_midpoint_tuple = [()] * (int((dim * (dim - 1)) / 2)) + curind = 0 + for ind1 in range(0, len(cur_index_to_node_tuple)): + for ind2 in range(ind1+1, len(cur_index_to_node_tuple)): + i = cur_index_to_node_tuple[ind1] + j = cur_index_to_node_tuple[ind2] + #print(i, j) + for k in range(0, dim-1): + cur = int((i[k] + j[k]) / 2) + cur_index_to_midpoint_tuple[curind] = cur_index_to_midpoint_tuple[curind] + (cur,) + curind += 1 + self.index_to_node_tuple.append(cur_index_to_node_tuple) + self.index_to_midpoint_tuple.append(cur_index_to_midpoint_tuple) + ''' + self.ray_vertices = np.empty([len(mesh.groups[0].vertex_indices) * + len(mesh.groups[0].vertex_indices[0]) * (len(mesh.groups[0].vertex_indices[0]) - 1) / 2, 2], + dtype=np.int32) + self.ray_elements = np.zeros([len(mesh.groups[0].vertex_indices) * + len(mesh.groups[0].vertex_indices[0]) * (len(mesh.groups[0].vertex_indices[0]) - 1) + / 2, 1, 2], dtype=np.int32) + self.vertex_to_ray = [] + for i in mesh.vertices[0]: + self.vertex_to_ray.append([]); + count = 0 + for grp in mesh.groups: + for i in range(0, len(grp.vertex_indices)): + for j in range(0, len(grp.vertex_indices[i])): + for k in range(j + 1, len(grp.vertex_indices[i])): + self.ray_vertices[count][0] = grp.vertex_indices[i][j] + self.ray_vertices[count][1] = grp.vertex_indices[i][k] + temp1 = VertexRay(count, 0) + self.vertex_to_ray[grp.vertex_indices[i][j]].append(temp1) + temp2 = VertexRay(count, 1) + self.vertex_to_ray[grp.vertex_indices[i][k]].append(temp2) + count += 1 + ind = 0 + #print(self.ray_vertices) + for i in self.ray_vertices: + which = 0 + for grp in mesh.groups: + for j in range(0, len(grp.vertex_indices)): + count = 0 + for k in grp.vertex_indices[j]: + if k == i[0] or k == i[1]: + count += 1 + #found an element sharing edge + if count == 2: + self.ray_elements[ind][0][which] = j + which += 1 + if which == 2: + break + ind += 1 + ''' + + # }}} + + # {{{ helper routines + + def get_refine_base_index(self): + if self.last_split_elements is None: + return 0 + else: + return self.last_mesh.nelements - len(self.last_split_elements) + + def get_empty_refine_flags(self): + return np.zeros( + 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 + #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): + queue = [cur_node] + res = [] + while queue: + vertex = queue.pop(0) + #if leaf node + if vertex.left is None and vertex.right is None: + res.append(vertex) + else: + queue.append(vertex.left) + queue.append(vertex.right) + return res + + def get_subtree(self, cur_node): + queue = [cur_node] + res = [] + while queue: + vertex = queue.pop(0) + res.append(vertex) + if not (vertex.left is None and vertex.right is None): + queue.append(vertex.left) + queue.append(vertex.right) + return res + + def apply_diff(self, cur_node, new_hanging_vertex_elements): + for el in cur_node.adjacent_add_diff: + if el not in cur_node.adjacent_elements: + cur_node.adjacent_elements.append(el) + if el not in new_hanging_vertex_elements[cur_node.left_vertex]: + new_hanging_vertex_elements[cur_node.left_vertex].append(el) + if el not in new_hanging_vertex_elements[cur_node.right_vertex]: + new_hanging_vertex_elements[cur_node.right_vertex].append(el) + if cur_node.left is not None and cur_node.right is not None: + if el not in cur_node.left.adjacent_add_diff: + cur_node.left.adjacent_add_diff.append(el) + if el not in cur_node.right.adjacent_add_diff: + cur_node.right.adjacent_add_diff.append(el) + cur_node.adjacent_add_diff = [] + +# def propagate(self, cur_node, new_hanging_vertex_elements): +# if cur_node.parent is not None: +# parent_node = cur_node.parent +# self.propagate(parent_node, new_hanging_vertex_elements) +# self.apply_diff(parent_node, new_hanging_vertex_elements) + + def get_root(self, cur_node): + while(cur_node.parent is not None): + cur_node = cur_node.parent + return cur_node + + def propagate_tree(self, cur_node, new_hanging_vertex_elements, element_to_element): + vertex_tuple = (min(cur_node.left_vertex, cur_node.right_vertex), max(cur_node.left_vertex, cur_node.right_vertex)) + self.seen_tuple[vertex_tuple] = True + self.apply_diff(cur_node, new_hanging_vertex_elements) + if cur_node.left is not None and cur_node.right is not None: + self.propagate_tree(cur_node.left, new_hanging_vertex_elements, element_to_element) + self.propagate_tree(cur_node.right, new_hanging_vertex_elements, element_to_element) + else: + for el in cur_node.adjacent_elements: + element_to_element[el].update(cur_node.adjacent_elements) + for el in new_hanging_vertex_elements[cur_node.left_vertex]: + element_to_element[el].update(new_hanging_vertex_elements[cur_node.left_vertex]) + for el in new_hanging_vertex_elements[cur_node.right_vertex]: + element_to_element[el].update(new_hanging_vertex_elements[cur_node.right_vertex]) + +# def remove_from_subtree(self, cur_node, new_hanging_vertex_elements, to_remove): +# if self.lazy: +# self.propagate(cur_node, new_hanging_vertex_elements) +# if to_remove in cur_node.adjacent_add_diff: +# cur_node.adjacent_add_diff.remove(to_remove) +# if to_remove not in cur_node.adjacent_remove_diff: +# cur_node.adjacent_remove_diff.append(to_remove) +# else: +# subtree = self.get_subtree(cur_node) +# for node in subtree: +# if to_remove in node.adjacent_elements: +# node.adjacent_elements.remove(to_remove) +# if to_remove in new_hanging_vertex_elements[node.left_vertex]: +# new_hanging_vertex_elements[node.left_vertex].remove(to_remove) +# if to_remove in new_hanging_vertex_elements[node.right_vertex]: +# new_hanging_vertex_elements[node.right_vertex].remove(to_remove) + + def add_to_subtree(self, cur_node, new_hanging_vertex_elements, to_add): + if self.lazy: + if to_add not in cur_node.adjacent_add_diff: + cur_node.adjacent_add_diff.append(to_add) + else: + subtree = self.get_subtree(cur_node) + for node in subtree: + if to_add not in node.adjacent_elements: + node.adjacent_elements.append(to_add) + if to_add not in new_hanging_vertex_elements[node.left_vertex]: + new_hanging_vertex_elements[node.left_vertex].append(to_add) + if to_add not in new_hanging_vertex_elements[node.right_vertex]: + new_hanging_vertex_elements[node.right_vertex].append(to_add) + + # }}} + + # {{{ refinement + + def refine(self, refine_flags): + """ + :arg refine_flags: a :class:`numpy.ndarray` of dtype bool of length ``mesh.nelements`` + indicating which elements should be split. + """ + + #vertices and groups for next generation + nvertices = len(self.last_mesh.vertices[0]) + + groups = [] + + midpoint_already = set() + grpn = 0 + totalnelements = 0 + + for grp in self.last_mesh.groups: + iel_base = grp.element_nr_base + nelements = 0 + for iel_grp in range(grp.nelements): + nelements += 1 + vertex_indices = grp.vertex_indices[iel_grp] + if refine_flags[iel_base+iel_grp]: + cur_dim = len(grp.vertex_indices[iel_grp])-1 + nelements += len(self.simplex_result[cur_dim]) - 1 + for i in range(len(vertex_indices)): + for j in range(i+1, len(vertex_indices)): + i_index = vertex_indices[i] + j_index = vertex_indices[j] + index_tuple = (i_index, j_index) if i_index < j_index else (j_index, i_index) + if index_tuple not in midpoint_already and \ + self.pair_map[index_tuple].midpoint is None: + nvertices += 1 + midpoint_already.add(index_tuple) + groups.append(np.empty([nelements, len(grp.vertex_indices[0])], dtype=np.int32)) + grpn += 1 + totalnelements += nelements + + vertices = np.empty([len(self.last_mesh.vertices), nvertices]) + + new_hanging_vertex_element = [ + [] for i in range(nvertices)] + +# def remove_element_from_connectivity(vertices, new_hanging_vertex_elements, to_remove): +# #print(vertices) +# import itertools +# if len(vertices) == 2: +# min_vertex = min(vertices[0], vertices[1]) +# max_vertex = max(vertices[0], vertices[1]) +# ray = self.pair_map[(min_vertex, max_vertex)] +# self.remove_from_subtree(ray, new_hanging_vertex_elements, to_remove) +# return +# +# cur_dim = len(vertices)-1 +# element_rays = [] +# midpoints = [] +# split_possible = True +# for i in range(len(vertices)): +# for j in range(i+1, len(vertices)): +# min_vertex = min(vertices[i], vertices[j]) +# max_vertex = max(vertices[i], vertices[j]) +# element_rays.append(self.pair_map[(min_vertex, max_vertex)]) +# if element_rays[len(element_rays)-1].midpoint is not None: +# midpoints.append(element_rays[len(element_rays)-1].midpoint) +# else: +# split_possible = False + + #for node in element_rays: + #self.remove_from_subtree(node, new_hanging_vertex_elements, to_remove) + #if split_possible: +# if split_possible: +# node_tuple_to_coord = {} +# for node_index, node_tuple in enumerate(self.index_to_node_tuple[cur_dim]): +# node_tuple_to_coord[node_tuple] = vertices[node_index] +# for midpoint_index, midpoint_tuple in enumerate(self.index_to_midpoint_tuple[cur_dim]): +# node_tuple_to_coord[midpoint_tuple] = midpoints[midpoint_index] +# for i in range(len(self.simplex_result[cur_dim])): +# next_vertices = [] +# for j in range(len(self.simplex_result[cur_dim][i])): +# next_vertices.append(node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]]) +# all_rays_present = True +# for v1 in range(len(next_vertices)): +# for v2 in range(v1+1, len(next_vertices)): +# vertex_tuple = (min(next_vertices[v1], next_vertices[v2]), max(next_vertices[v1], next_vertices[v2])) +# if vertex_tuple not in self.pair_map: +# all_rays_present = False +# if all_rays_present: +# remove_element_from_connectivity(next_vertices, new_hanging_vertex_elements, to_remove) +# else: +# split_possible = False +# if not split_possible: +# next_vertices_list = list(itertools.combinations(vertices, len(vertices)-1)) +# 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]) + max_vertex = max(vertices[0], vertices[1]) + ray = self.pair_map[(min_vertex, max_vertex)] + self.add_to_subtree(ray, new_hanging_vertex_elements, to_add) + return + + cur_dim = len(vertices)-1 + element_rays = [] + midpoints = [] + split_possible = True + for i in range(len(vertices)): + for j in range(i+1, len(vertices)): + min_vertex = min(vertices[i], vertices[j]) + max_vertex = max(vertices[i], vertices[j]) + element_rays.append(self.pair_map[(min_vertex, max_vertex)]) + if element_rays[len(element_rays)-1].midpoint is not None: + midpoints.append(element_rays[len(element_rays)-1].midpoint) + else: + split_possible = False + #for node in element_rays: + #self.add_to_subtree(node, new_hanging_vertex_elements, to_add) + if split_possible: + node_tuple_to_coord = {} + for node_index, node_tuple in enumerate(self.index_to_node_tuple[cur_dim]): + node_tuple_to_coord[node_tuple] = vertices[node_index] + for midpoint_index, midpoint_tuple in enumerate(self.index_to_midpoint_tuple[cur_dim]): + node_tuple_to_coord[midpoint_tuple] = midpoints[midpoint_index] + for i in range(len(self.simplex_result[cur_dim])): + next_vertices = [] + for j in range(len(self.simplex_result[cur_dim][i])): + next_vertices.append(node_tuple_to_coord[self.simplex_node_tuples[cur_dim][self.simplex_result[cur_dim][i][j]]]) + all_rays_present = True + for v1 in range(len(next_vertices)): + for v2 in range(v1+1, len(next_vertices)): + vertex_tuple = (min(next_vertices[v1], next_vertices[v2]), max(next_vertices[v1], next_vertices[v2])) + if vertex_tuple not in self.pair_map: + all_rays_present = False + if all_rays_present: + add_element_to_connectivity(next_vertices, new_hanging_vertex_elements, to_add) + else: + split_possible = False + if not split_possible: + next_vertices_list = list(itertools.combinations(vertices, len(vertices)-1)) + for next_vertices in next_vertices_list: + add_element_to_connectivity(next_vertices, new_hanging_vertex_elements, to_add) +# for node in element_rays: +# self.add_element_to_connectivity(node, new_hanging_vertex_elements, to_add) + # leaves = self.get_subtree(node) + # for leaf in leaves: + # if to_add not in leaf.adjacent_elements: + # leaf.adjacent_elements.append(to_add) + # if to_add not in new_hanging_vertex_elements[leaf.left_vertex]: + # new_hanging_vertex_elements[leaf.left_vertex].append(to_add) + # if to_add not in new_hanging_vertex_elements[leaf.right_vertex]: + # new_hanging_vertex_elements[leaf.right_vertex].append(to_add) + +# next_element_rays = [] +# for i in range(len(element_rays)): +# for j in range(i+1, len(element_rays)): +# if element_rays[i].midpoint is not None and element_rays[j].midpoint is not None: +# min_midpoint = min(element_rays[i].midpoint, element_rays[j].midpoint) +# max_midpoint = max(element_rays[i].midpoint, element_rays[j].midpoint) +# vertex_pair = (min_midpoint, max_midpoint) +# if vertex_pair in self.pair_map: +# next_element_rays.append(self.pair_map[vertex_pair]) +# cur_next_rays = [] +# if element_rays[i].left_vertex == element_rays[j].left_vertex: +# cur_next_rays = [element_rays[i].left, element_rays[j].left, self.pair_map[vertex_pair]] +# if element_rays[i].right_vertex == element_rays[j].right_vertex: +# cur_next_rays = [element_rays[i].right, element_rays[j].right, self.pair_map[vertex_pair]] +# if element_rays[i].left_vertex == element_rays[j].right_vertex: +# cur_next_rays = [element_rays[i].left, element_rays[j].right, self.pair_map[vertex_pair]] +# if element_rays[i].right_vertex == element_rays[j].left_vertex: +# cur_next_rays = [element_rays[i].right, element_rays[j].left, self.pair_map[vertex_pair]] +# assert (cur_next_rays != []) +# #print cur_next_rays +# add_element_to_connectivity(cur_next_rays, new_hanging_vertex_elements, to_add) +# else: +# return +# else: +# 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 + for iel_grp in range(nelements_in_grp): + vertex_indices = grp[iel_grp] + 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]) + max_index = max(vertex_indices[i], vertex_indices[j]) + cur_node = self.pair_map[(min_index, max_index)] + #print iel_base+iel_grp, cur_node.left_vertex, cur_node.right_vertex + #if (iel_base + iel_grp) not in cur_node.adjacent_elements: + #print min_index, max_index + #print iel_base + iel_grp, cur_node.left_vertex, cur_node.right_vertex, cur_node.adjacent_elements + #assert (4 in new_hanging_vertex_elements[cur_node.left_vertex] or 4 in new_hanging_vertex_elements[cur_node.right_vertex]) + assert ((iel_base+iel_grp) in cur_node.adjacent_elements) + 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] + import copy + if i == 0: + new_hanging_vertex_element[j] = copy.deepcopy(self.hanging_vertex_element[j]) + grpn = 0 + for grp in self.last_mesh.groups: + for iel_grp in range(grp.nelements): + for i in range(len(grp.vertex_indices[iel_grp])): + groups[grpn][iel_grp][i] = grp.vertex_indices[iel_grp][i] + grpn += 1 + + grpn = 0 + vertices_index = len(self.last_mesh.vertices[0]) + nelements_in_grp = grp.nelements + 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]) + max_index = max(vertex_indices[i], vertex_indices[j]) + cur_node = self.pair_map[(min_index, max_index)] + if cur_node.midpoint is None: + cur_node.midpoint = vertices_index + import copy + cur_node.left = TreeRayNode(min_index, vertices_index, + copy.deepcopy(cur_node.adjacent_elements)) + cur_node.left.parent = cur_node + cur_node.right = TreeRayNode(max_index, vertices_index, + copy.deepcopy(cur_node.adjacent_elements)) + cur_node.right.parent = cur_node + vertex_pair1 = (min_index, vertices_index) + 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_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)): + for j in range(i+1, len(midpoint_vertices)): + min_index = min(midpoint_vertices[i], midpoint_vertices[j]) + max_index = max(midpoint_vertices[i], midpoint_vertices[j]) + vertex_pair = (min_index, max_index) + if vertex_pair in self.pair_map: + continue + self.pair_map[vertex_pair] = TreeRayNode(min_index, max_index, []) + 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] + 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])): + 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: + #quadrilateral +# 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] +# 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: + iel_base = grp.element_nr_base + for iel_grp in range(grp.nelements): + for i in range(len(grp.vertex_indices[iel_grp])): + for j in range(i+1, len(grp.vertex_indices[iel_grp])): + min_vert = min(grp.vertex_indices[iel_grp][i], grp.vertex_indices[iel_grp][j]) + max_vert = max(grp.vertex_indices[iel_grp][i], grp.vertex_indices[iel_grp][j]) + vertex_pair = (min_vert, max_vert) + root_ray = self.get_root(self.pair_map[vertex_pair]) + if root_ray not in self.seen_tuple: + self.seen_tuple[root_ray] = True + cur_tree = self.get_subtree(root_ray) + for node in cur_tree: + node.adjacent_elements = [] + new_hanging_vertex_element[node.left_vertex] = [] + new_hanging_vertex_element[node.right_vertex] = [] + + self.seen_tuple.clear() + + nelements_in_grp = grp.nelements + for grp in groups: + for iel_grp in range(len(grp)): + add_verts = [] + for i in range(len(grp[iel_grp])): + add_verts.append(grp[iel_grp][i]) + add_element_to_connectivity(add_verts, new_hanging_vertex_element, iel_base+iel_grp) + #assert ray connectivity + #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)) + + from meshmode.mesh import Mesh + + self.previous_mesh = self.last_mesh + self.last_mesh = Mesh( + vertices, grp, + nodal_adjacency=self.generate_nodal_adjacency( + totalnelements, nvertices, groups), + vertex_id_dtype=self.last_mesh.vertex_id_dtype, + element_id_dtype=self.last_mesh.element_id_dtype) + return self.last_mesh + + # }}} + + def print_rays(self, ind): + for i in range(len(self.last_mesh.groups[0].vertex_indices[ind])): + for j in range(i+1, len(self.last_mesh.groups[0].vertex_indices[ind])): + mn = min(self.last_mesh.groups[0].vertex_indices[ind][i], + self.last_mesh.groups[0].vertex_indices[ind][j]) + mx = max(self.last_mesh.groups[0].vertex_indices[ind][j], + self.last_mesh.groups[0].vertex_indices[ind][i]) + vertex_pair = (mn, mx) + print('LEFT VERTEX:', self.pair_map[vertex_pair].left_vertex) + print('RIGHT VERTEX:', self.pair_map[vertex_pair].right_vertex) + print('ADJACENT:') + rays = self.get_leaves(self.pair_map[vertex_pair]) + for k in rays: + print(k.adjacent_elements) + ''' + def print_rays(self, groups, ind): + for i in range(len(groups[0][ind])): + for j in range(i+1, len(groups[0][ind])): + mn = min(groups[0][ind][i], groups[0][ind][j]) + mx = max(groups[0][ind][i], groups[0][ind][j]) + vertex_pair = (mn, mx) + print('LEFT VERTEX:', self.pair_map[vertex_pair].left_vertex) + print('RIGHT VERTEX:', self.pair_map[vertex_pair].right_vertex) + print('ADJACENT:') + rays = self.get_leaves(self.pair_map[vertex_pair]) + for k in rays: + print(k.adjacent_elements) + ''' + + def print_hanging_elements(self, ind): + for i in self.last_mesh.groups[0].vertex_indices[ind]: + print("IND:", i, self.hanging_vertex_element[i]) + + # {{{ generate adjacency + + def generate_nodal_adjacency(self, nelements, nvertices, groups): + # medium-term FIXME: make this an incremental update + # rather than build-from-scratch + vertex_to_element = [[] for i in range(nvertices)] + element_index = 0 + for grp in groups: + for iel_grp in range(len(grp)): + for ivertex in grp[iel_grp]: + vertex_to_element[ivertex].append(element_index) + element_index += 1 + element_to_element = [set() for i in range(nelements)] + element_index = 0 + if self.lazy: + for grp in groups: + for iel_grp in range(len(grp)): + for i in range(len(grp[iel_grp])): + for j in range(i+1, len(grp[iel_grp])): + vertex_pair = (min(grp[iel_grp][i], grp[iel_grp][j]), max(grp[iel_grp][i], grp[iel_grp][j])) + #print 'iel:', iel_grp, 'pair:', vertex_pair + if vertex_pair not in self.seen_tuple: + self.propagate_tree(self.get_root(self.pair_map[vertex_pair]), self.hanging_vertex_element, element_to_element) + #print self.pair_map[vertex_pair].left_vertex, self.pair_map[vertex_pair].right_vertex, self.pair_map[vertex_pair].adjacent_elements, self.hanging_vertex_element[self.pair_map[vertex_pair].left_vertex], self.hanging_vertex_element[self.pair_map[vertex_pair].right_vertex] + + else: + for grp in groups: + for iel_grp in range(len(grp)): + for ivertex in grp[iel_grp]: + element_to_element[element_index].update( + vertex_to_element[ivertex]) + if self.hanging_vertex_element[ivertex]: + for hanging_element in self.hanging_vertex_element[ivertex]: + if element_index != hanging_element: + element_to_element[element_index].update([hanging_element]) + element_to_element[hanging_element].update([element_index]) + for i in range(len(grp[iel_grp])): + for j in range(i+1, len(grp[iel_grp])): + vertex_pair = (min(grp[iel_grp][i], grp[iel_grp][j]), + max(grp[iel_grp][i], grp[iel_grp][j])) + #element_to_element[element_index].update( + #self.pair_map[vertex_pair].adjacent_elements) + queue = [self.pair_map[vertex_pair]] + while queue: + vertex = queue.pop(0) + #if leaf node + if vertex.left is None and vertex.right is None: + assert(element_index in vertex.adjacent_elements) + element_to_element[element_index].update( + vertex.adjacent_elements) + else: + queue.append(vertex.left) + queue.append(vertex.right) + ''' + if self.hanging_vertex_element[ivertex] and element_index != self.hanging_vertex_element[ivertex][0]: + element_to_element[element_index].update([self.hanging_vertex_element[ivertex][0]]) + element_to_element[self.hanging_vertex_element[ivertex][0]].update([element_index]) + ''' + element_index += 1 + print(len(element_to_element)) + for iel, neighbors in enumerate(element_to_element): + if iel in neighbors: + neighbors.remove(iel) + #print(self.ray_elements) + ''' + for ray in self.rays: + curnode = ray.first + while curnode is not None: + if len(curnode.value.elements) >= 2: + if curnode.value.elements[0] is not None: + element_to_element[curnode.value.elements[0]].update(curnode.value.elements) + if curnode.value.elements[1] is not None: + element_to_element[curnode.value.elements[1]].update(curnode.value.elements) + if len(curnode.value.velements) >= 2: + if curnode.value.velements[0] is not None: + element_to_element[curnode.value.velements[0]].update(curnode.value.velements) + if curnode.value.velements[1] is not None: + element_to_element[curnode.value.velements[1]].update(curnode.value.velements) + curnode = curnode.next + ''' + ''' + for i in self.ray_elements: + for j in i: + #print j[0], j[1] + element_to_element[j[0]].update(j) + element_to_element[j[1]].update(j) + ''' + #print element_to_element + lengths = [len(el_list) for el_list in element_to_element] + neighbors_starts = np.cumsum( + np.array([0] + lengths, dtype=self.last_mesh.element_id_dtype), + # cumsum silently widens integer types + dtype=self.last_mesh.element_id_dtype) + from pytools import flatten + neighbors = np.array( + list(flatten(element_to_element)), + dtype=self.last_mesh.element_id_dtype) + + assert neighbors_starts[-1] == len(neighbors) + + from meshmode.mesh import NodalAdjacency + return NodalAdjacency(neighbor_starts=neighbors_starts, neighbors=neighbors) + + # }}} + + +# vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py new file mode 100644 index 0000000000000000000000000000000000000000..9a78c31d6cd2dbe5bdb2f3729815b06630c96c7f --- /dev/null +++ b/meshmode/mesh/refinement/utils.py @@ -0,0 +1,121 @@ +from __future__ import division, absolute_import, print_function + +__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.neighbor_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 = np.linalg.solve(transformation, vertex_transformed) + + is_connected = ( + 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/meshmode/mesh/tesselate.py b/meshmode/mesh/tesselate.py new file mode 100644 index 0000000000000000000000000000000000000000..5d0e03aef4b9b2a8ac9ead1571dbdea1290ea2a7 --- /dev/null +++ b/meshmode/mesh/tesselate.py @@ -0,0 +1,74 @@ +from pytools import generate_nonnegative_integer_tuples_summing_to_at_most \ + as gnitstam + + +def add_tuples(a, b): + return tuple(ac+bc for ac, bc in zip(a, b)) + + +def tesselatetri(): + result = [] + + node_tuples = list(gnitstam(2, 2)) + node_dict = 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 cube into six tets. + # subtets 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 = 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] diff --git a/setup.py b/setup.py index 0b15c94320e2073c270f68b247422c0cbbb1b844..06d1d63ced0aad542ba53440b716073828ab8eb5 100644 --- a/setup.py +++ b/setup.py @@ -42,6 +42,7 @@ def main(): packages=find_packages(), install_requires=[ "numpy", + "llist", "modepy", "meshpy>=2014.1", "six", diff --git a/test/test_meshmode.py b/test/test_meshmode.py index b3ca2022c55989281a9637339cd5d0cee94a81dc..01e6c35cd95cd5d87270d43962d57bd25b605514 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -751,7 +751,7 @@ def test_as_python(): # }}} -# {{{ lookup tree for element finding +# {{{ test lookup tree for element finding def test_lookup_tree(do_plot=False): from meshmode.mesh.generation import make_curve_mesh, cloverleaf diff --git a/test/test_refinement.py b/test/test_refinement.py new file mode 100644 index 0000000000000000000000000000000000000000..9c36b9750ebe87a899a8dd5c5f17f4191ed1b0b4 --- /dev/null +++ b/test/test_refinement.py @@ -0,0 +1,221 @@ +from __future__ import division, absolute_import, print_function + +__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 pytest +import pyopencl as cl +import pyopencl.clmath # noqa + +import numpy as np +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl + as pytest_generate_tests) +from meshmode.mesh.generation import ( # noqa + generate_icosahedron, generate_box_mesh) +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(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;" % h] + ) + + +def random_refine_flags(fract, mesh): + all_els = list(range(mesh.nelements)) + + flags = np.zeros(mesh.nelements) + from random import shuffle + shuffle(all_els) + for i in range(int(mesh.nelements * fract)): + flags[all_els[i]] = 1 + + return flags + + +def uniform_refine_flags(mesh): + return np.ones(mesh.nelements) + + +@pytest.mark.parametrize(("case_name", "mesh_gen", "flag_gen", "num_generations"), [ + # Fails? + # ("icosahedron", + # partial(generate_icosahedron, 1, order=1), + # partial(random_refine_flags, 0.4), + # 3), + + ("rect2d_rand", + partial(generate_box_mesh, ( + np.linspace(0, 1, 3), + np.linspace(0, 1, 3), + ), order=1), + partial(random_refine_flags, 0.4), + 4), + + ("rect2d_unif", + partial(generate_box_mesh, ( + np.linspace(0, 1, 2), + np.linspace(0, 1, 2), + ), order=1), + uniform_refine_flags, + 3), + + ("blob2d_rand", + gen_blob_mesh, + partial(random_refine_flags, 0.4), + 4), + + ("rect3d_rand", + partial(generate_box_mesh, ( + np.linspace(0, 1, 2), + np.linspace(0, 1, 3), + np.linspace(0, 1, 2), + ), order=1), + partial(random_refine_flags, 0.4), + 3), + + ("rect3d_unif", + partial(generate_box_mesh, ( + np.linspace(0, 1, 2), + np.linspace(0, 1, 2)), order=1), + uniform_refine_flags, + 3), + ]) +def test_refinement(case_name, mesh_gen, flag_gen, num_generations): + from random import seed + seed(13) + + mesh = mesh_gen() + + r = Refiner(mesh) + + for igen in range(num_generations): + flags = flag_gen(mesh) + mesh = r.refine(flags) + + 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: + exec(sys.argv[1]) + else: + from py.test.cmdline import main + main([__file__]) + +# vim: fdm=marker