diff --git a/meshmode/discretization/connection/__init__.py b/meshmode/discretization/connection/__init__.py index 30e886eda68462ba6ecd0222e74c247f6279050c..bd4f1292a083ade8273f5f92e01efebd66d652e8 100644 --- a/meshmode/discretization/connection/__init__.py +++ b/meshmode/discretization/connection/__init__.py @@ -35,12 +35,14 @@ from meshmode.discretization.connection.same_mesh import \ make_same_mesh_connection from meshmode.discretization.connection.face import ( FACE_RESTR_INTERIOR, FACE_RESTR_ALL, - make_face_restriction, make_face_to_all_faces_embedding) + 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 - +from meshmode.discretization.connection.chained import \ + flatten_chained_connection import logging logger = logging.getLogger(__name__) @@ -53,11 +55,13 @@ __all__ = [ "make_face_restriction", "make_face_to_all_faces_embedding", "make_opposite_face_connection", - "make_refinement_connection" + "make_refinement_connection", + "flatten_chained_connection" ] __doc__ = """ .. autoclass:: DiscretizationConnection +.. autoclass:: ChainedDiscretizationConnection .. autoclass:: DirectDiscretizationConnection .. autofunction:: make_same_mesh_connection @@ -72,6 +76,8 @@ __doc__ = """ .. autofunction:: make_refinement_connection +.. autofunction:: flatten_chained_connection + Implementation details ^^^^^^^^^^^^^^^^^^^^^^ diff --git a/meshmode/discretization/connection/chained.py b/meshmode/discretization/connection/chained.py new file mode 100644 index 0000000000000000000000000000000000000000..221ce2d37650f3110dc82e9981b7d44d5a74569b --- /dev/null +++ b/meshmode/discretization/connection/chained.py @@ -0,0 +1,260 @@ +from __future__ import division, print_function, absolute_import + +__copyright__ = """Copyright (C) 2018 Alexandru Fikl""" + +__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 + +from pytools import Record + +import modepy as mp + + +# {{{ flatten chained connection + +class _ConnectionBatchData(Record): + pass + + +def _iterbatches(groups): + for igrp, grp in enumerate(groups): + for ibatch, batch in enumerate(grp.batches): + yield (igrp, ibatch), (grp, batch) + + +def _build_element_lookup_table(queue, conn): + el_table = [np.full(g.nelements, -1, dtype=np.int) + for g in conn.to_discr.groups] + + for (igrp, _), (_, batch) in _iterbatches(conn.groups): + el_table[igrp][batch.to_element_indices.get(queue)] = \ + batch.from_element_indices.get(queue) + + return el_table + + +def _build_new_group_table(from_conn, to_conn): + def find_batch(nodes, gtb): + for igrp, batches in enumerate(gtb): + for ibatch, batch in enumerate(batches): + if np.allclose(nodes, batch.result_unit_nodes): + return (igrp, ibatch) + return (-1, -1) + + nfrom_groups = len(from_conn.groups) + nto_groups = len(to_conn.groups) + + # construct a table from (old groups) -> (new groups) + # NOTE: we try to reduce the number of new groups and batches by matching + # the `result_unit_nodes` and only adding a new batch if necessary + grp_to_grp = {} + batch_info = [[] for i in range(nfrom_groups * nto_groups)] + for (igrp, ibatch), (fgrp, fbatch) in _iterbatches(from_conn.groups): + for (jgrp, jbatch), (tgrp, tbatch) in _iterbatches(to_conn.groups): + # compute result_unit_nodes + ffgrp = from_conn.from_discr.groups[fbatch.from_group_index] + from_matrix = mp.resampling_matrix( + ffgrp.basis(), + fbatch.result_unit_nodes, + ffgrp.unit_nodes) + result_unit_nodes = from_matrix.dot(ffgrp.unit_nodes.T) + + tfgrp = to_conn.from_discr.groups[tbatch.from_group_index] + to_matrix = mp.resampling_matrix( + tfgrp.basis(), + tbatch.result_unit_nodes, + tfgrp.unit_nodes) + result_unit_nodes = to_matrix.dot(result_unit_nodes).T + + # find new (group, batch) + (igrp_new, ibatch_new) = find_batch(result_unit_nodes, batch_info) + if igrp_new < 0: + igrp_new = nto_groups * igrp + jgrp + ibatch_new = len(batch_info[igrp_new]) + + batch_info[igrp_new].append(_ConnectionBatchData( + from_group_index=fbatch.from_group_index, + result_unit_nodes=result_unit_nodes, + to_element_face=tbatch.to_element_face)) + + grp_to_grp[igrp, ibatch, jgrp, jbatch] = (igrp_new, ibatch_new) + + return grp_to_grp, batch_info + + +def _build_batches(queue, from_bins, to_bins, batch): + from meshmode.discretization.connection import \ + InterpolationBatch + + def to_device(x): + return cl.array.to_device(queue, np.asarray(x)) + + for ibatch, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins)): + yield InterpolationBatch( + from_group_index=batch[ibatch].from_group_index, + from_element_indices=to_device(from_bin), + to_element_indices=to_device(to_bin), + result_unit_nodes=batch[ibatch].result_unit_nodes, + to_element_face=batch[ibatch].to_element_face) + + +def flatten_chained_connection(queue, connection): + """Collapse a connection into a direct connection. + + If the given connection is already a + :class:`~meshmode.discretization.connection.DirectDiscretizationConnection` + nothing is done. However, if the connection is a + :class:`~meshmode.discretization.connection.ChainedDiscretizationConnection`, + a new direct connection is constructed that transports from + :attr:`connection.from_discr` to :attr:`connection.to_discr`. + + The new direct connection will have a number of groups and batches that + is, at worse, the product of all the connections in the chain. For + example, if we consider a connection between a discretization and a + two-level refinement, both levels will have :math:`n` groups and + :math:`m + 1` batches per group, where :math:`m` is the number of + subdivisions of an element (exact number depends on implementation + details in + :func:`~meshmode.discretizationc.connection.make_refinement_connection`). + However, a direct connection from level :math:`0` to level :math:`2` + will have at worst :math:`n^2` groups and each group will have + :math:`(m + 1)^2` batches. + + .. warning:: + + If a large number of connections is chained, the number of groups and + batches can become very large. + + :arg queue: An instance of :class:`pyopencl.CommandQueue`. + :arg connection: An instance of + :class:`~meshmode.discretization.connection.DiscretizationConnection`. + :return: An instance of + :class:`~meshmode.discretization.connection.DirectDiscretizationConnection`. + """ + from meshmode.discretization.connection import ( + DirectDiscretizationConnection, + DiscretizationConnectionElementGroup, + make_same_mesh_connection) + + if not hasattr(connection, 'connections'): + return connection + + if not connection.connections: + return make_same_mesh_connection(connection.to_discr, + connection.from_discr) + + # recursively build direct connections + connections = connection.connections + direct_connections = [] + for conn in connections: + direct_connections.append(flatten_chained_connection(queue, conn)) + + # merge all the direct connections + from_conn = direct_connections[0] + for to_conn in direct_connections[1:]: + el_table = _build_element_lookup_table(queue, from_conn) + grp_to_grp, batch_info = _build_new_group_table(from_conn, to_conn) + + # distribute the indices to new groups and batches + from_bins = [[np.empty(0, dtype=np.int) for _ in g] for g in batch_info] + to_bins = [[np.empty(0, dtype=np.int) for _ in g] for g in batch_info] + + for (igrp, ibatch), (_, from_batch) in _iterbatches(from_conn.groups): + from_to_element_indices = from_batch.to_element_indices.get(queue) + + for (jgrp, jbatch), (_, to_batch) in _iterbatches(to_conn.groups): + igrp_new, ibatch_new = grp_to_grp[igrp, ibatch, jgrp, jbatch] + + jfrom = to_batch.from_element_indices.get(queue) + jto = to_batch.to_element_indices.get(queue) + + mask = np.isin(jfrom, from_to_element_indices) + from_bins[igrp_new][ibatch_new] = \ + np.hstack([from_bins[igrp_new][ibatch_new], + el_table[igrp][jfrom[mask]]]) + to_bins[igrp_new][ibatch_new] = \ + np.hstack([to_bins[igrp_new][ibatch_new], + jto[mask]]) + + # build new groups + groups = [] + for igrp, (from_bin, to_bin) in enumerate(zip(from_bins, to_bins)): + groups.append(DiscretizationConnectionElementGroup( + list(_build_batches(queue, from_bin, to_bin, + batch_info[igrp])))) + + from_conn = DirectDiscretizationConnection( + from_discr=from_conn.from_discr, + to_discr=to_conn.to_discr, + groups=groups, + is_surjective=connection.is_surjective) + + return from_conn + +# }}} + + +# {{{ build chained resample matrix + +def make_full_resample_matrix(queue, connection): + """Build a dense matrix representing the discretization connection. + + This is based on + :func:`~meshmode.discretization.connection.DirectDiscretizationConnection.full_resample_matrix`. + If a chained connection is given, the matrix is constructed recursively + for each connection and multiplied left to right. + + .. warning:: + + This method will be very slow, both in terms of speed and memory + usage, and should only be used for testing or if absolutely necessary. + + :arg queue: a :class:`pyopencl.CommandQueue`. + :arg connection: a + :class:`~meshmode.discretization.connection.DiscretizationConnection`. + :return: a :class:`pyopencl.array.Array` of shape + `(connection.from_discr.nnodes, connection.to_discr.nnodes)`. + """ + + if hasattr(connection, "full_resample_matrix"): + return connection.full_resample_matrix(queue) + + if not hasattr(connection, 'connections'): + raise TypeError('connection is not chained') + + if not connection.connections: + result = np.eye(connection.to_discr.nnodes) + return cl.array.to_device(queue, result) + + acc = make_full_resample_matrix(queue, connection.connections[0]).get(queue) + for conn in connection.connections[1:]: + resampler = make_full_resample_matrix(queue, conn).get(queue) + acc = resampler.dot(acc) + + return cl.array.to_device(queue, acc) + +# }}} + + +# vim: foldmethod=marker diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py index 4413117d0ea69aef016c353b24a16011cc55c9be..7ab02e701413df407f8d3c9f2711ec855639d99e 100644 --- a/meshmode/mesh/refinement/no_adjacency.py +++ b/meshmode/mesh/refinement/no_adjacency.py @@ -134,7 +134,7 @@ class RefinerWithoutAdjacency(object): # }}} def refine_uniformly(self): - flags = np.ones(self._last_mesh.nelements, dtype=bool) + flags = np.ones(self._current_mesh.nelements, dtype=bool) self.refine(flags) # {{{ refinement top-level diff --git a/test/test_chained.py b/test/test_chained.py new file mode 100644 index 0000000000000000000000000000000000000000..c0f5a09e320ce4442f1747a35c06718f4067d1dd --- /dev/null +++ b/test/test_chained.py @@ -0,0 +1,343 @@ +from __future__ import division, absolute_import, print_function + +__copyright__ = "Copyright (C) 2018 Alexandru Fikl" + +__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 pyopencl.clmath # noqa + +import pytest +from pyopencl.tools import ( # noqa + pytest_generate_tests_for_pyopencl + as pytest_generate_tests) + +import logging +logger = logging.getLogger(__name__) + + +def create_discretization(queue, ndim, + nelements=42, + mesh_order=5, + discr_order=5): + ctx = queue.context + + # construct mesh + if ndim == 2: + from functools import partial + from meshmode.mesh.generation import make_curve_mesh, ellipse + mesh = make_curve_mesh(partial(ellipse, 2), + np.linspace(0.0, 1.0, nelements + 1), + order=mesh_order) + elif ndim == 3: + from meshmode.mesh.generation import generate_torus + mesh = generate_torus(10.0, 5.0, order=mesh_order) + else: + raise ValueError("Unsupported dimension: {}".format(ndim)) + + # create discretization + from meshmode.discretization import Discretization + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + discr = Discretization(ctx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(discr_order)) + + return discr + + +def create_refined_connection(queue, discr, threshold=0.3): + from meshmode.mesh.refinement import RefinerWithoutAdjacency + from meshmode.discretization.connection import make_refinement_connection + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + flags = np.random.rand(discr.mesh.nelements) < threshold + refiner = RefinerWithoutAdjacency(discr.mesh) + refiner.refine(flags) + + discr_order = discr.groups[0].order + connection = make_refinement_connection(refiner, discr, + InterpolatoryQuadratureSimplexGroupFactory(discr_order)) + + return connection + + +def create_face_connection(queue, discr): + from meshmode.discretization.connection import FACE_RESTR_ALL + from meshmode.discretization.connection import make_face_restriction + from meshmode.discretization.poly_element import \ + InterpolatoryQuadratureSimplexGroupFactory + + discr_order = discr.groups[0].order + connection = make_face_restriction(discr, + InterpolatoryQuadratureSimplexGroupFactory(discr_order), + FACE_RESTR_ALL, + per_face_groups=True) + + return connection + + +@pytest.mark.skip(reason='implementation detail') +@pytest.mark.parametrize("ndim", [2, 3]) +def test_chained_batch_table(ctx_factory, ndim, visualize=False): + from meshmode.discretization.connection.chained import \ + _build_element_lookup_table + + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + discr = create_discretization(queue, ndim) + connections = [] + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + + from meshmode.discretization.connection import ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) + + conn = chained.connections[0] + el_table = _build_element_lookup_table(queue, conn) + for igrp, grp in enumerate(conn.groups): + for ibatch, batch in enumerate(grp.batches): + ifrom = batch.from_element_indices.get(queue) + jfrom = el_table[igrp][batch.to_element_indices.get(queue)] + + assert np.all(ifrom == jfrom) + assert np.min(el_table[igrp]) >= 0 + + +@pytest.mark.skip(reason='implementation detail') +@pytest.mark.parametrize("ndim", [2, 3]) +def test_chained_new_group_table(ctx_factory, ndim, visualize=False): + from meshmode.discretization.connection.chained import \ + _build_new_group_table + + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + discr = create_discretization(queue, ndim, + nelements=8, + mesh_order=2, + discr_order=2) + connections = [] + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + + grp_to_grp, grp_info = _build_new_group_table(connections[0], + connections[1]) + if visualize: + import matplotlib.pyplot as pt + + pt.figure(figsize=(10, 8)) + for k, v in grp_to_grp.items(): + print(k) + print(v) + + igrp, ibatch, jgrp, jbatch = k + mgroup, mbatch = v + from_group_index = connections[0].groups[igrp] \ + .batches[ibatch].from_group_index + + from_unit_nodes = connections[0].from_discr \ + .groups[from_group_index].unit_nodes + to_unit_nodes = grp_info[mgroup][mbatch].result_unit_nodes + + if ndim == 2: + pt.plot(from_unit_nodes, 'o') + pt.plot(to_unit_nodes, '^') + else: + pt.plot(from_unit_nodes[0], from_unit_nodes[1], 'o') + pt.plot(to_unit_nodes[0], to_unit_nodes[1], '^') + + pt.savefig('test_grp_to_grp_{}d_{:05d}_{:05d}.png' + .format(ndim, mgroup, mbatch), dpi=300) + pt.clf() + + +@pytest.mark.parametrize("ndim", [2, 3]) +def test_chained_connection(ctx_factory, ndim, visualize=False): + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + discr = create_discretization(queue, ndim, + nelements=10, + mesh_order=5, + discr_order=5) + connections = [] + conn = create_refined_connection(queue, discr, threshold=np.inf) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr, threshold=np.inf) + connections.append(conn) + + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) + + def f(x): + from six.moves import reduce + return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + + x = connections[0].from_discr.nodes().with_queue(queue) + fx = f(x) + f1 = chained(queue, fx).get(queue) + f2 = connections[1](queue, connections[0](queue, fx)).get(queue) + + assert np.allclose(f1, f2) + + +@pytest.mark.skip(reason='slow test') +@pytest.mark.parametrize("ndim", [2, 3]) +def test_chained_full_resample_matrix(ctx_factory, ndim, visualize=False): + from meshmode.discretization.connection.chained import \ + make_full_resample_matrix + + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + discr = create_discretization(queue, ndim, + mesh_order=2, + discr_order=2) + connections = [] + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) + + def f(x): + from six.moves import reduce + return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + + resample_mat = make_full_resample_matrix(queue, chained).get(queue) + + x = connections[0].from_discr.nodes().with_queue(queue) + fx = f(x) + f1 = np.dot(resample_mat, fx.get(queue)) + f2 = chained(queue, fx).get(queue) + f3 = connections[1](queue, connections[0](queue, fx)).get(queue) + + assert np.allclose(f1, f2) + assert np.allclose(f2, f3) + + +@pytest.mark.parametrize(("ndim", "chain_type"), [ + (2, 1), (2, 2), (3, 1), (3, 3)]) +def test_chained_to_direct(ctx_factory, ndim, chain_type, + nelements=128, visualize=False): + import time + from meshmode.discretization.connection.chained import \ + flatten_chained_connection + + ctx = ctx_factory() + queue = cl.CommandQueue(ctx) + + discr = create_discretization(queue, ndim, nelements=nelements) + connections = [] + if chain_type == 1: + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + elif chain_type == 2: + conn = create_refined_connection(queue, discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + conn = create_refined_connection(queue, conn.to_discr) + connections.append(conn) + elif chain_type == 3 and ndim == 3: + conn = create_refined_connection(queue, discr, threshold=np.inf) + connections.append(conn) + conn = create_face_connection(queue, conn.to_discr) + connections.append(conn) + else: + raise ValueError('unknown test case') + + from meshmode.discretization.connection import \ + ChainedDiscretizationConnection + chained = ChainedDiscretizationConnection(connections) + + t_start = time.time() + direct = flatten_chained_connection(queue, chained) + t_end = time.time() + if visualize: + print('[TIME] Flatten: {:.5e}'.format(t_end - t_start)) + + if chain_type < 3: + to_element_indices = np.full(direct.to_discr.mesh.nelements, 0, + dtype=np.int) + for grp in direct.groups: + for batch in grp.batches: + for i in batch.to_element_indices.get(queue): + to_element_indices[i] += 1 + assert np.min(to_element_indices) > 0 + + def f(x): + from six.moves import reduce + return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) + + x = connections[0].from_discr.nodes().with_queue(queue) + fx = f(x) + + t_start = time.time() + f1 = direct(queue, fx).get(queue) + t_end = time.time() + if visualize: + print('[TIME] Direct: {:.5e}'.format(t_end - t_start)) + + t_start = time.time() + f2 = chained(queue, fx).get(queue) + t_end = time.time() + if visualize: + print('[TIME] Chained: {:.5e}'.format(t_end - t_start)) + + if visualize and ndim == 2: + import matplotlib.pyplot as pt + + pt.figure(figsize=(10, 8)) + pt.plot(f1, label='Direct') + pt.plot(f2, label='Chained') + pt.ylim([np.min(f2) - 0.1, np.max(f2) + 0.1]) + pt.legend() + pt.savefig('test_chained_to_direct.png', dpi=300) + pt.clf() + + assert np.allclose(f1, f2) + + +if __name__ == "__main__": + import sys + if len(sys.argv) > 1: + exec(sys.argv[1]) + else: + from pytest import main + main([__file__]) + +# vim: fdm=marker diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 7173fc29d359b9a6e7ef65f889f819e3290a455a..29dc2ed1075aa31621e9d4d5ac1578c9732f8bd8 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -1039,65 +1039,6 @@ def test_quad_multi_element(): plt.show() -# {{{ ChainedDiscretizationConnection - -def test_ChainedDiscretizationConnection(ctx_getter): # noqa - mesh_order = 5 - order = 5 - npanels = 10 - group_factory = InterpolatoryQuadratureSimplexGroupFactory - - def refine_flags(mesh): - return np.ones(mesh.nelements) - - cl_ctx = ctx_getter() - queue = cl.CommandQueue(cl_ctx) - - from functools import partial - from meshmode.discretization import Discretization - from meshmode.discretization.connection import make_refinement_connection - from meshmode.mesh.generation import make_curve_mesh, ellipse - - mesh = make_curve_mesh( - partial(ellipse, 1), np.linspace(0, 1, npanels + 1), - order=mesh_order) - - discr = Discretization(cl_ctx, mesh, group_factory(order)) - - connections = [] - - from meshmode.mesh.refinement import Refiner - refiner = Refiner(mesh) - - def refine_discr(discr): - mesh = discr.mesh - flags = refine_flags(mesh) - refiner.refine(flags) - connections.append( - make_refinement_connection(refiner, discr, group_factory(order))) - return connections[-1].to_discr - - discr = refine_discr(discr) - refine_discr(discr) - - from meshmode.discretization.connection import ( - ChainedDiscretizationConnection) - - chained_conn = ChainedDiscretizationConnection(connections) - - def f(x): - from six.moves import reduce - return 0.1 * reduce(lambda x, y: x * cl.clmath.sin(5 * y), x) - - x = connections[0].from_discr.nodes().with_queue(queue) - - assert np.allclose( - chained_conn(queue, f(x)).get(queue), - connections[1](queue, connections[0](queue, f(x))).get(queue)) - -# }}} - - if __name__ == "__main__": import sys if len(sys.argv) > 1: