From 62604f23721e55aa5def2141f4e8475e24921016 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 9 Jun 2020 23:30:16 -0500 Subject: [PATCH 01/13] add visualization support for VTK_LAGRANGE elements --- meshmode/discretization/visualization.py | 447 +++++++++++++++-------- 1 file changed, 292 insertions(+), 155 deletions(-) diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index e985f3ff..d0d5b965 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -37,7 +37,7 @@ __doc__ = """ """ -# {{{ visualizer +# {{{ helpers def separate_by_real_and_imag(data, real_only): for name, field in data: @@ -70,6 +70,19 @@ def separate_by_real_and_imag(data, real_only): yield (name, field) +def resample_and_get(queue, conn, vec): + from pytools.obj_array import with_object_array_or_scalar + + def resample_and_get_one(fld): + from numbers import Number + if isinstance(fld, Number): + return np.ones(conn.to_discr.nnodes) * fld + else: + return conn(queue, fld).get(queue=queue) + + return with_object_array_or_scalar(resample_and_get_one, vec) + + class _VisConnectivityGroup(Record): """ .. attribute:: vis_connectivity @@ -97,117 +110,128 @@ class _VisConnectivityGroup(Record): def primitive_element_size(self): return self.vis_connectivity.shape[2] +# }}} -class Visualizer(object): + +# {{{ vtk visualizers + +class VTKVisualizer(object): """ - .. automethod:: show_scalar_in_mayavi - .. automethod:: show_scalar_in_matplotlib_3d .. automethod:: write_vtk_file """ def __init__(self, connection, element_shrink_factor=None): + if element_shrink_factor is None: + element_shrink_factor = 1.0 + + self.element_shrink_factor = element_shrink_factor self.connection = connection self.discr = connection.from_discr self.vis_discr = connection.to_discr - if element_shrink_factor is None: - element_shrink_factor = 1 + @property + def version(self): + return "0.1" - self.element_shrink_factor = element_shrink_factor + # {{{ connectivity - def _resample_and_get(self, queue, vec): - from pytools.obj_array import with_object_array_or_scalar + @property + def simplex_cell_types(self): + import pyvisfile.vtk as vtk + return { + 1: vtk.VTK_LINE, + 2: vtk.VTK_TRIANGLE, + 3: vtk.VTK_TETRA, + } - def resample_and_get_one(fld): - from numbers import Number - if isinstance(fld, Number): - return np.ones(self.connection.to_discr.nnodes) * fld - else: - return self.connection(queue, fld).get(queue=queue) + @property + def tensor_cell_types(self): + import pyvisfile.vtk as vtk + return { + 1: vtk.VTK_LINE, + 2: vtk.VTK_QUAD, + 3: vtk.VTK_HEXAHEDRON, + } + + def connectivity_for_element_group(self, grp): + from pytools import ( + generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, + generate_nonnegative_integer_tuples_below as gnitb) + from meshmode.mesh import TensorProductElementGroup, SimplexElementGroup + + if isinstance(grp.mesh_el_group, SimplexElementGroup): + node_tuples = list(gnitstam(grp.order, grp.dim)) + + from modepy.tools import submesh + el_connectivity = np.array( + submesh(node_tuples), + dtype=np.intp) + + vtk_cell_type = self.simplex_cell_types[grp.dim] + + elif isinstance(grp.mesh_el_group, TensorProductElementGroup): + node_tuples = list(gnitb(grp.order+1, grp.dim)) + node_tuple_to_index = dict( + (nt, i) for i, nt in enumerate(node_tuples)) + + def add_tuple(a, b): + return tuple(ai+bi for ai, bi in zip(a, b)) + + el_offsets = { + 1: [(0,), (1,)], + 2: [(0, 0), (1, 0), (1, 1), (0, 1)], + 3: [ + (0, 0, 0), + (1, 0, 0), + (1, 1, 0), + (0, 1, 0), + (0, 0, 1), + (1, 0, 1), + (1, 1, 1), + (0, 1, 1), + ] + }[grp.dim] + + el_connectivity = np.array([ + [ + node_tuple_to_index[add_tuple(origin, offset)] + for offset in el_offsets] + for origin in gnitb(grp.order, grp.dim)]) + + vtk_cell_type = self.tensor_cell_types[grp.dim] - return with_object_array_or_scalar(resample_and_get_one, vec) + else: + raise NotImplementedError("visualization for element groups " + "of type '%s'" % type(grp.mesh_el_group).__name__) - # {{{ vis sub-element connectivity + assert len(node_tuples) == grp.nunit_nodes + return el_connectivity, vtk_cell_type + @property @memoize_method - def _vis_connectivity(self): + def cells(self): + return np.hstack([ + vgrp.vis_connectivity.reshape(-1) for vgrp in self.groups + ]) + + @property + @memoize_method + def groups(self): """ :return: a list of :class:`_VisConnectivityGroup` instances. """ # Assume that we're using modepy's default node ordering. - from pytools import ( - generate_nonnegative_integer_tuples_summing_to_at_most as gnitstam, - generate_nonnegative_integer_tuples_below as gnitb) - from meshmode.mesh import TensorProductElementGroup, SimplexElementGroup - result = [] - - from pyvisfile.vtk import ( - VTK_LINE, VTK_TRIANGLE, VTK_TETRA, - VTK_QUAD, VTK_HEXAHEDRON) - subel_nr_base = 0 - for group in self.vis_discr.groups: - if isinstance(group.mesh_el_group, SimplexElementGroup): - node_tuples = list(gnitstam(group.order, group.dim)) - - from modepy.tools import submesh - el_connectivity = np.array( - submesh(node_tuples), - dtype=np.intp) - vtk_cell_type = { - 1: VTK_LINE, - 2: VTK_TRIANGLE, - 3: VTK_TETRA, - }[group.dim] - - elif isinstance(group.mesh_el_group, TensorProductElementGroup): - node_tuples = list(gnitb(group.order+1, group.dim)) - node_tuple_to_index = dict( - (nt, i) for i, nt in enumerate(node_tuples)) - - def add_tuple(a, b): - return tuple(ai+bi for ai, bi in zip(a, b)) - - el_offsets = { - 1: [(0,), (1,)], - 2: [(0, 0), (1, 0), (1, 1), (0, 1)], - 3: [ - (0, 0, 0), - (1, 0, 0), - (1, 1, 0), - (0, 1, 0), - (0, 0, 1), - (1, 0, 1), - (1, 1, 1), - (0, 1, 1), - ] - }[group.dim] - - el_connectivity = np.array([ - [ - node_tuple_to_index[add_tuple(origin, offset)] - for offset in el_offsets] - for origin in gnitb(group.order, group.dim)]) - - vtk_cell_type = { - 1: VTK_LINE, - 2: VTK_QUAD, - 3: VTK_HEXAHEDRON, - }[group.dim] + for grp in self.vis_discr.groups: + el_connectivity, vtk_cell_type = \ + self.connectivity_for_element_group(grp) - else: - raise NotImplementedError("visualization for element groups " - "of type '%s'" % type(group.mesh_el_group).__name__) - - assert len(node_tuples) == group.nunit_nodes - vis_connectivity = ( - group.node_nr_base + np.arange( - 0, group.nelements*group.nunit_nodes, group.nunit_nodes - )[:, np.newaxis, np.newaxis] - + el_connectivity).astype(np.intp) + offsets = grp.node_nr_base \ + + np.arange(0, grp.nnodes, grp.nunit_nodes).reshape(-1, 1, 1) + vis_connectivity = (offsets + el_connectivity).astype(np.intp) vgrp = _VisConnectivityGroup( vis_connectivity=vis_connectivity, @@ -221,6 +245,181 @@ class Visualizer(object): # }}} + def write_vtk_file(self, file_name, names_and_fields, + compressor=None, real_only=False, overwrite=False): + from pyvisfile.vtk import ( + UnstructuredGrid, DataArray, + AppendedDataXMLGenerator, + VF_LIST_OF_COMPONENTS) + + with cl.CommandQueue(self.vis_discr.cl_context) as queue: + nodes = self.vis_discr.nodes().get(queue) + + names_and_fields = [ + (name, resample_and_get(queue, self.connection, f)) + for name, f in names_and_fields + ] + + # {{{ create cell_types + + nsubelements = sum(vgrp.nsubelements for vgrp in self.groups) + cell_types = np.empty(nsubelements, dtype=np.uint8) + cell_types.fill(255) + + for vgrp in self.groups: + isubelements = np.s_[ + vgrp.subelement_nr_base: + vgrp.subelement_nr_base + vgrp.nsubelements] + cell_types[isubelements] = vgrp.vtk_cell_type + + assert (cell_types < 255).all() + + # }}} + + # {{{ shrink elements + + if abs(self.element_shrink_factor - 1.0) > 1.0e-14: + for vgrp in self.vis_discr.groups: + nodes_view = vgrp.view(nodes) + el_centers = np.mean(nodes_view, axis=-1) + nodes_view[:] = ( + (self.element_shrink_factor * nodes_view) + + (1-self.element_shrink_factor) + * el_centers[:, :, np.newaxis]) + + # }}} + + # {{{ create grid + + points = DataArray("points", + nodes.reshape(self.vis_discr.ambient_dim, -1), + vector_format=VF_LIST_OF_COMPONENTS) + + grid = UnstructuredGrid( + (self.vis_discr.nnodes, points), + cells=self.cells, + cell_types=cell_types) + + for name, field in separate_by_real_and_imag(names_and_fields, real_only): + grid.add_pointdata( + DataArray(name, field, vector_format=VF_LIST_OF_COMPONENTS) + ) + + # }}} + + # {{{ write + + import os + from meshmode import FileExistsError + if os.path.exists(file_name): + if overwrite: + os.remove(file_name) + else: + raise FileExistsError("output file '%s' already exists" % file_name) + + with open(file_name, "w") as outf: + generator = AppendedDataXMLGenerator( + compressor=compressor, + vtk_file_version=self.version) + + generator(grid).write(outf) + + # }}} + + +class VTKLagrangeVisualizer(VTKVisualizer): + @property + def version(self): + # NOTE: version 2.2 has an updated ordering for the hexahedron + # elements that is not supported currently + # https://gitlab.kitware.com/vtk/vtk/-/merge_requests/6678 + return "2.1" + + @property + def simplex_cell_types(self): + import pyvisfile.vtk as vtk + return { + 1: vtk.VTK_LAGRANGE_CURVE, + 2: vtk.VTK_LAGRANGE_TRIANGLE, + 3: vtk.VTK_LAGRANGE_TETRAHEDRON, + } + + @property + def tensor_cell_types(self): + import pyvisfile.vtk as vtk + return { + 1: vtk.VTK_LAGRANGE_CURVE, + 2: vtk.VTK_LAGRANGE_QUADRILATERAL, + 3: vtk.VTK_LAGRANGE_HEXAHEDRON, + } + + def connectivity_for_element_group(self, grp): + from meshmode.mesh import SimplexElementGroup + + if isinstance(grp.mesh_el_group, SimplexElementGroup): + from pyvisfile.vtk.vtk_ordering import ( + vtk_lagrange_simplex_node_tuples, + vtk_lagrange_simplex_node_tuples_to_permutation) + + node_tuples = vtk_lagrange_simplex_node_tuples(grp.dim, grp.order) + el_connectivity = np.array( + vtk_lagrange_simplex_node_tuples_to_permutation(node_tuples), + dtype=np.intp).reshape(1, 1, -1) + + vtk_cell_type = self.simplex_cell_types[grp.dim] + + else: + raise NotImplementedError("visualization for element groups " + "of type '%s'" % type(grp.mesh_el_group).__name__) + + assert len(node_tuples) == grp.nunit_nodes + return el_connectivity, vtk_cell_type + + @property + @memoize_method + def cells(self): + connectivity = np.hstack([ + grp.vis_connectivity.reshape(-1) + for grp in self.groups + ]) + offsets = np.hstack([ + np.arange(grp.nunit_nodes, grp.nnodes + 1, grp.nunit_nodes) + for grp in self.vis_discr.groups + ]) + + from pyvisfile.vtk import DataArray + return ( + self.vis_discr.mesh.nelements, + DataArray("connectivity", connectivity), + DataArray("offsets", offsets) + ) + +# }}} + + +# {{{ visualizer + +class Visualizer(object): + """ + .. automethod:: show_scalar_in_mayavi + .. automethod:: show_scalar_in_matplotlib_3d + .. automethod:: write_vtk_file + """ + + def __init__(self, connection, + element_shrink_factor=None, + use_high_order_vtk=False): + self.connection = connection + self.discr = connection.from_discr + self.vis_discr = connection.to_discr + + if use_high_order_vtk: + self.vtk = VTKLagrangeVisualizer(connection, + element_shrink_factor=element_shrink_factor) + else: + self.vtk = VTKVisualizer(connection, + element_shrink_factor=element_shrink_factor) + # {{{ mayavi def show_scalar_in_mayavi(self, field, **kwargs): @@ -231,7 +430,7 @@ class Visualizer(object): with cl.CommandQueue(self.vis_discr.cl_context) as queue: nodes = self.vis_discr.nodes().with_queue(queue).get() - field = self._resample_and_get(queue, field) + field = resample_and_get(queue, self.connection, field) assert nodes.shape[0] == self.vis_discr.ambient_dim #mlab.points3d(nodes[0], nodes[1], 0*nodes[0]) @@ -280,74 +479,12 @@ class Visualizer(object): compressor=None, real_only=False, overwrite=False): + self.vtk.write_vtk_file(file_name, names_and_fields, + compressor=compressor, + real_only=real_only, + overwrite=overwrite) - from pyvisfile.vtk import ( - UnstructuredGrid, DataArray, - AppendedDataXMLGenerator, - VF_LIST_OF_COMPONENTS) - - with cl.CommandQueue(self.vis_discr.cl_context) as queue: - nodes = self.vis_discr.nodes().with_queue(queue).get() - - names_and_fields = [ - (name, self._resample_and_get(queue, fld)) - for name, fld in names_and_fields] - - vc_groups = self._vis_connectivity() - - # {{{ create cell_types - - nsubelements = sum(vgrp.nsubelements for vgrp in vc_groups) - cell_types = np.empty(nsubelements, dtype=np.uint8) - cell_types.fill(255) - for vgrp in vc_groups: - cell_types[ - vgrp.subelement_nr_base: - vgrp.subelement_nr_base + vgrp.nsubelements] = \ - vgrp.vtk_cell_type - assert (cell_types < 255).all() - - # }}} - - if self.element_shrink_factor != 1: - for vgrp in self.vis_discr.groups: - nodes_view = vgrp.view(nodes) - el_centers = np.mean(nodes_view, axis=-1) - nodes_view[:] = ( - (self.element_shrink_factor * nodes_view) - + (1-self.element_shrink_factor) - * el_centers[:, :, np.newaxis]) - - grid = UnstructuredGrid( - (self.vis_discr.nnodes, - DataArray("points", - nodes.reshape(self.vis_discr.ambient_dim, -1), - vector_format=VF_LIST_OF_COMPONENTS)), - cells=np.hstack([ - vgrp.vis_connectivity.reshape(-1) - for vgrp in vc_groups]), - cell_types=cell_types) - - # for name, field in separate_by_real_and_imag(cell_data, real_only): - # grid.add_celldata(DataArray(name, field, - # vector_format=VF_LIST_OF_COMPONENTS)) - - for name, field in separate_by_real_and_imag(names_and_fields, real_only): - grid.add_pointdata(DataArray(name, field, - vector_format=VF_LIST_OF_COMPONENTS)) - - import os - from meshmode import FileExistsError - if os.path.exists(file_name): - if overwrite: - os.remove(file_name) - else: - raise FileExistsError("output file '%s' already exists" % file_name) - - with open(file_name, "w") as outf: - AppendedDataXMLGenerator(compressor)(grid).write(outf) - - # }}} + # }}} # {{{ matplotlib 3D @@ -365,7 +502,7 @@ class Visualizer(object): with cl.CommandQueue(self.vis_discr.cl_context) as queue: nodes = self.vis_discr.nodes().with_queue(queue).get() - field = self._resample_and_get(queue, field) + field = resample_and_get(queue, self.connection, field) assert nodes.shape[0] == self.vis_discr.ambient_dim -- GitLab From 4358b053c42ff288ffd2c187d20c1545ef9d27c5 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 10 Jun 2020 11:32:05 -0500 Subject: [PATCH 02/13] fix some visualizer bugs --- meshmode/discretization/poly_element.py | 20 ++++++--- meshmode/discretization/visualization.py | 32 +++++++++----- test/test_meshmode.py | 56 +++++++++++------------- 3 files changed, 59 insertions(+), 49 deletions(-) diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index 436b1fb1..c9063b20 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -42,6 +42,7 @@ Group types .. autoclass:: PolynomialWarpAndBlendElementGroup .. autoclass:: PolynomialEquidistantSimplexElementGroup .. autoclass:: LegendreGaussLobattoTensorProductElementGroup +.. autoclass:: EquidistantTensorProductElementGroup Group factories ^^^^^^^^^^^^^^^ @@ -276,16 +277,21 @@ class LegendreGaussLobattoTensorProductElementGroup(PolynomialElementGroupBase): @memoize_method def from_mesh_interp_matrix(self): - from modepy.modes import tensor_product_basis, jacobi - from functools import partial meg = self.mesh_el_group + return mp.resampling_matrix( + self.basis(), + self.unit_nodes, + meg.unit_nodes) - basis = tensor_product_basis( - self.dim, tuple( - partial(jacobi, 0, 0, i) - for i in range(meg.order+1))) - return mp.resampling_matrix(basis, self.unit_nodes, meg.unit_nodes) +class EquidistantTensorProductElementGroup( + LegendreGaussLobattoTensorProductElementGroup): + @property + @memoize_method + def unit_nodes(self): + from modepy.nodes import tensor_product_nodes, equidistant_nodes + return tensor_product_nodes( + self.dim, equidistant_nodes(1, self.order)) # }}} diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index d0d5b965..b5c9d97d 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -333,7 +333,7 @@ class VTKLagrangeVisualizer(VTKVisualizer): # NOTE: version 2.2 has an updated ordering for the hexahedron # elements that is not supported currently # https://gitlab.kitware.com/vtk/vtk/-/merge_requests/6678 - return "2.1" + return "2.0" @property def simplex_cell_types(self): @@ -383,7 +383,8 @@ class VTKLagrangeVisualizer(VTKVisualizer): for grp in self.groups ]) offsets = np.hstack([ - np.arange(grp.nunit_nodes, grp.nnodes + 1, grp.nunit_nodes) + grp.node_nr_base + + np.arange(grp.nunit_nodes, grp.nnodes + 1, grp.nunit_nodes) for grp in self.vis_discr.groups ]) @@ -557,26 +558,35 @@ class Visualizer(object): # }}} -def make_visualizer(queue, discr, vis_order, element_shrink_factor=None): +def make_visualizer(queue, discr, vis_order, + element_shrink_factor=None, use_high_order_vtk=False): from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import ( - PolynomialWarpAndBlendElementGroup, - LegendreGaussLobattoTensorProductElementGroup, - OrderAndTypeBasedGroupFactory) + + if use_high_order_vtk: + from meshmode.discretization.poly_element import ( + PolynomialEquidistantSimplexElementGroup as SimplexElementGroup, + EquidistantTensorProductElementGroup as TensorElementGroup) + else: + from meshmode.discretization.poly_element import ( + PolynomialWarpAndBlendElementGroup as SimplexElementGroup, + LegendreGaussLobattoTensorProductElementGroup as TensorElementGroup) + + from meshmode.discretization.poly_element import OrderAndTypeBasedGroupFactory vis_discr = Discretization( discr.cl_context, discr.mesh, OrderAndTypeBasedGroupFactory( vis_order, - simplex_group_class=PolynomialWarpAndBlendElementGroup, - tensor_product_group_class=( - LegendreGaussLobattoTensorProductElementGroup)), + simplex_group_class=SimplexElementGroup, + tensor_product_group_class=TensorElementGroup), real_dtype=discr.real_dtype) + from meshmode.discretization.connection import \ make_same_mesh_connection return Visualizer( make_same_mesh_connection(vis_discr, discr), - element_shrink_factor=element_shrink_factor) + element_shrink_factor=element_shrink_factor, + use_high_order_vtk=use_high_order_vtk) # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 047362d2..8e740c77 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -548,7 +548,7 @@ def test_element_orientation(): ("ball", lambda: mgen.generate_icosahedron(1, 1)), ("torus", lambda: mgen.generate_torus(5, 1)), ]) -def test_3d_orientation(ctx_factory, what, mesh_gen_func, visualize=False): +def test_orientation_3d(ctx_factory, what, mesh_gen_func, visualize=False): pytest.importorskip("pytential") logging.basicConfig(level=logging.INFO) @@ -562,7 +562,7 @@ def test_3d_orientation(ctx_factory, what, mesh_gen_func, visualize=False): from meshmode.discretization import Discretization discr = Discretization(ctx, mesh, - PolynomialWarpAndBlendGroupFactory(1)) + PolynomialWarpAndBlendGroupFactory(3)) from pytential import bind, sym @@ -592,9 +592,9 @@ def test_3d_orientation(ctx_factory, what, mesh_gen_func, visualize=False): if visualize: from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, discr, 1) + vis = make_visualizer(queue, discr, 3) - vis.write_vtk_file("normals.vtu", [ + vis.write_vtk_file("orientation_3d_normals_%s.vtu" % what, [ ("normals", normals), ]) @@ -634,7 +634,7 @@ def test_merge_and_map(ctx_factory, visualize=False): from meshmode.mesh.processing import merge_disjoint_meshes, affine_map mesh2 = affine_map(mesh, A=np.eye(mesh.ambient_dim), - b=np.array([5, 0, 0])[:mesh.ambient_dim]) + b=np.array([2, 0, 0])[:mesh.ambient_dim]) mesh3 = merge_disjoint_meshes((mesh2, mesh)) mesh3.facial_adjacency_groups @@ -650,7 +650,7 @@ def test_merge_and_map(ctx_factory, visualize=False): from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(queue, discr, 3, element_shrink_factor=0.8) - vis.write_vtk_file("merged.vtu", []) + vis.write_vtk_file("merge_and_map.vtu", []) # }}} @@ -715,20 +715,15 @@ def test_sanity_single_element(ctx_factory, dim, order, visualize=False): # }}} - # {{{ visualizers - - from meshmode.discretization.visualization import make_visualizer - #vol_vis = make_visualizer(queue, vol_discr, 4) - bdry_vis = make_visualizer(queue, bdry_discr, 4) - - # }}} - from pytential import bind, sym bdry_normals = bind(bdry_discr, sym.normal(dim))(queue).as_vector(dtype=object) if visualize: - bdry_vis.write_vtk_file("boundary.vtu", [ - ("bdry_normals", bdry_normals) + from meshmode.discretization.visualization import make_visualizer + bdry_vis = make_visualizer(queue, bdry_discr, 4) + + bdry_vis.write_vtk_file("sanity_single_element_boundary.vtu", [ + ("normals", bdry_normals) ]) normal_outward_check = bind(bdry_discr, @@ -799,8 +794,7 @@ def test_sanity_qhull_nd(ctx_factory, dim, order): ("ball-radius-1.step", 3), ]) @pytest.mark.parametrize("mesh_order", [1, 2]) -def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, - visualize=False): +def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, visualize=False): pytest.importorskip("pytential") logging.basicConfig(level=logging.INFO) @@ -841,14 +835,6 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, # }}} - # {{{ visualizers - - from meshmode.discretization.visualization import make_visualizer - vol_vis = make_visualizer(queue, vol_discr, 20) - bdry_vis = make_visualizer(queue, bdry_discr, 20) - - # }}} - from math import gamma true_surf = 2*np.pi**(dim/2)/gamma(dim/2) true_vol = true_surf/dim @@ -879,14 +865,22 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, print("SURF", true_surf, comp_surf) if visualize: - vol_vis.write_vtk_file("volume-h=%g.vtu" % h, [ + from meshmode.discretization.visualization import make_visualizer + vol_vis = make_visualizer(queue, vol_discr, 7) + bdry_vis = make_visualizer(queue, bdry_discr, 7) + + name = src_file.split("-")[0] + vol_vis.write_vtk_file("sanity_balls_volume_%s_%g.vtu" % (name, h), [ ("f", vol_one), ("area_el", bind( vol_discr, sym.area_element(mesh.ambient_dim, mesh.ambient_dim)) (queue)), ]) - bdry_vis.write_vtk_file("boundary-h=%g.vtu" % h, [("f", bdry_one)]) + + bdry_vis.write_vtk_file("sanity_balls_boundary_%s_%g.vtu" % (name, h), [ + ("f", bdry_one) + ]) # {{{ check normals point outward @@ -942,7 +936,7 @@ def test_box_mesh(ctx_factory, visualize=False): from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(queue, discr, 1) - vis.write_vtk_file("box.vtu", []) + vis.write_vtk_file("box_mesh.vtu", []) # }}} @@ -1142,7 +1136,7 @@ def test_vtk_overwrite(ctx_getter): import os from meshmode import FileExistsError - filename = "test_vtk_overwrite.vtu" + filename = "vtk_overwrite_temp.vtu" if os.path.exists(filename): os.remove(filename) @@ -1395,7 +1389,7 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(queue, discr, vis_order=order) - vis.write_vtk_file("test_mesh_multiple_groups.vtu", [ + vis.write_vtk_file("mesh_multiple_groups.vtu", [ ("group_id", group_id) ], overwrite=True) -- GitLab From d1ce9f41cff6a692736831fdc1cfa7a2eb54623a Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Wed, 17 Jun 2020 18:37:58 -0500 Subject: [PATCH 03/13] pull pyvisfile from git for examples --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 3c9da790..3ad15daa 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -57,7 +57,7 @@ Python 3 POCL Examples: - export PY_EXE=python3 - export PYOPENCL_TEST=portable:pthread # cython is here because pytential (for now, for TS) depends on it - - export EXTRA_INSTALL="pybind11 cython numpy mako pyvisfile matplotlib" + - export EXTRA_INSTALL="pybind11 cython numpy mako git+git://github.com/inducer/pyvisfile matplotlib" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project-and-run-examples.sh - ". ./build-py-project-and-run-examples.sh" tags: -- GitLab From 455959f8124941708c7d9b2c0bed0f2a80d759c2 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sun, 28 Jun 2020 22:26:07 -0500 Subject: [PATCH 04/13] fix renaming merge issues --- meshmode/discretization/visualization.py | 13 ++++++++----- test/test_meshmode.py | 18 +++++++++--------- 2 files changed, 17 insertions(+), 14 deletions(-) diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index 84056635..032b3418 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -77,7 +77,7 @@ def resample_to_numpy(conn, vec): and vec.dtype.char == "O" and not isinstance(vec, DOFArray)): from pytools.obj_array import obj_array_vectorize - return obj_array_vectorize(resample_to_numpy, vec) + return obj_array_vectorize(lambda x: resample_to_numpy(conn, x), vec) from numbers import Number if isinstance(vec, Number): @@ -168,9 +168,9 @@ class VTKVisualizer(object): if isinstance(grp.mesh_el_group, SimplexElementGroup): node_tuples = list(gnitstam(grp.order, grp.dim)) - from modepy.tools import submesh + from modepy.tools import simplex_submesh el_connectivity = np.array( - submesh(node_tuples), + simplex_submesh(node_tuples), dtype=np.intp) vtk_cell_type = self.simplex_cell_types[grp.dim] @@ -296,7 +296,10 @@ class VTKVisualizer(object): if abs(self.element_shrink_factor - 1.0) > 1.0e-14: node_nr_base = 0 for vgrp in self.vis_discr.groups: - nodes_view = nodes[:, node_nr_base:node_nr_base + vgrp.ndofs] + nodes_view = ( + nodes[:, node_nr_base:node_nr_base + vgrp.ndofs] + .reshape(nodes.shape[0], vgrp.nelements, vgrp.nunit_dofs)) + el_centers = np.mean(nodes_view, axis=-1) nodes_view[:] = ( (self.element_shrink_factor * nodes_view) @@ -410,7 +413,7 @@ class VTKLagrangeVisualizer(VTKVisualizer): grp.nunit_dofs, grp.nelements * grp.nunit_dofs + 1, grp.nunit_dofs) - for grp_offset, grp in enumerate(grp_offsets, self.vis_discr.groups) + for grp_offset, grp in zip(grp_offsets, self.vis_discr.groups) ]) from pyvisfile.vtk import DataArray diff --git a/test/test_meshmode.py b/test/test_meshmode.py index daadc29d..9b1ec987 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -654,11 +654,11 @@ def test_merge_and_map(ctx_factory, visualize=False): from meshmode.discretization import Discretization cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) - - discr = Discretization(cl_ctx, mesh3, discr_grp_factory) + actx = PyOpenCLArrayContext(queue) + discr = Discretization(actx, mesh3, discr_grp_factory) from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, discr, 3, element_shrink_factor=0.8) + vis = make_visualizer(actx, discr, 3, element_shrink_factor=0.8) vis.write_vtk_file("merge_and_map.vtu", []) # }}} @@ -943,12 +943,13 @@ def test_box_mesh(ctx_factory, visualize=False): PolynomialWarpAndBlendGroupFactory cl_ctx = ctx_factory() queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) - discr = Discretization(cl_ctx, mesh, - PolynomialWarpAndBlendGroupFactory(1)) + discr = Discretization(actx, mesh, + PolynomialWarpAndBlendGroupFactory(7)) from meshmode.discretization.visualization import make_visualizer - vis = make_visualizer(queue, discr, 1) + vis = make_visualizer(actx, discr, 7) vis.write_vtk_file("box_mesh.vtu", []) # }}} @@ -1401,9 +1402,8 @@ def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False): if visualize: group_id = discr.empty(actx, dtype=np.int) - for igrp, grp in enumerate(discr.groups): - group_id_view = grp.view(group_id) - group_id_view.fill(igrp) + for igrp, vec in enumerate(group_id): + vec.fill(igrp) from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(actx, discr, vis_order=order) -- GitLab From 261bf777799b3fa63c4e64f318f6b18160700fff Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Sat, 4 Jul 2020 15:12:01 -0500 Subject: [PATCH 05/13] update ci --- .gitlab-ci.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 3ad15daa..3c9da790 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -57,7 +57,7 @@ Python 3 POCL Examples: - export PY_EXE=python3 - export PYOPENCL_TEST=portable:pthread # cython is here because pytential (for now, for TS) depends on it - - export EXTRA_INSTALL="pybind11 cython numpy mako git+git://github.com/inducer/pyvisfile matplotlib" + - export EXTRA_INSTALL="pybind11 cython numpy mako pyvisfile matplotlib" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-py-project-and-run-examples.sh - ". ./build-py-project-and-run-examples.sh" tags: -- GitLab From e466d7f2814d4d4469c5aa53de750baee384bc99 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Tue, 14 Jul 2020 19:44:33 -0500 Subject: [PATCH 06/13] use visualize instead of do_plot --- test/test_meshmode.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 5a77567f..bfb6b539 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -54,9 +54,9 @@ logger = logging.getLogger(__name__) # {{{ circle mesh -def test_circle_mesh(do_plot=False): +def test_circle_mesh(visualize=False): from meshmode.mesh.io import generate_gmsh, FileSource - print("BEGIN GEN") + logger.info("BEGIN GEN") mesh = generate_gmsh( FileSource("circle.step"), 2, order=2, force_ambient_dim=2, @@ -64,18 +64,22 @@ def test_circle_mesh(do_plot=False): "-string", "Mesh.CharacteristicLengthMax = 0.05;"], target_unit="MM", ) - print("END GEN") - print(mesh.nelements) + logger.info("END GEN") + logger.info("nelements: %d", mesh.nelements) from meshmode.mesh.processing import affine_map mesh = affine_map(mesh, A=3*np.eye(2)) - if do_plot: + if visualize: from meshmode.mesh.visualization import draw_2d_mesh - draw_2d_mesh(mesh, fill=None, draw_nodal_adjacency=True, + draw_2d_mesh(mesh, + fill=None, + draw_vertex_numbers=False, + draw_nodal_adjacency=True, set_bounding_box=True) import matplotlib.pyplot as pt - pt.show() + pt.axis("equal") + pt.savefig("circle_mesh", dpi=300) # }}} @@ -926,11 +930,11 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, visualize=False): # {{{ rect/box mesh generation -def test_rect_mesh(do_plot=False): +def test_rect_mesh(visualize=False): from meshmode.mesh.generation import generate_regular_rect_mesh mesh = generate_regular_rect_mesh() - if do_plot: + if visualize: from meshmode.mesh.visualization import draw_2d_mesh draw_2d_mesh(mesh, fill=None, draw_nodal_adjacency=True) import matplotlib.pyplot as pt @@ -991,7 +995,7 @@ def test_as_python(): # {{{ test lookup tree for element finding -def test_lookup_tree(do_plot=False): +def test_lookup_tree(visualize=False): from meshmode.mesh.generation import make_curve_mesh, cloverleaf mesh = make_curve_mesh(cloverleaf, np.linspace(0, 1, 1000), order=3) @@ -1009,7 +1013,7 @@ def test_lookup_tree(do_plot=False): for igrp, iel in tree.generate_matches(pt): print(igrp, iel) - if do_plot: + if visualize: with open("tree.dat", "w") as outf: tree.visualize(outf) -- GitLab From f474394e5cb57a7c228563ef85bdd1e88afc2897 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 17 Jul 2020 21:01:02 -0500 Subject: [PATCH 07/13] rename VTKVisualizer to VTKConnectivity --- meshmode/discretization/visualization.py | 276 +++++++++++------------ test/test_meshmode.py | 55 +++++ 2 files changed, 191 insertions(+), 140 deletions(-) diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index 032b3418..19f4f44c 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -93,7 +93,7 @@ class _VisConnectivityGroup(Record): """ .. attribute:: vis_connectivity - an array of shape ``(group.nelements,nsubelements,primitive_element_size)`` + an array of shape ``(group.nelements, nsubelements, primitive_element_size)`` .. attribute:: vtk_cell_type @@ -121,16 +121,15 @@ class _VisConnectivityGroup(Record): # {{{ vtk visualizers -class VTKVisualizer(object): - """ - .. automethod:: write_vtk_file - """ +class VTKConnectivity: + """Connectivity for standard linear VTK element types. - def __init__(self, connection, element_shrink_factor=None): - if element_shrink_factor is None: - element_shrink_factor = 1.0 + .. attribute:: version + .. attribute:: cells + .. attribute:: groups + """ - self.element_shrink_factor = element_shrink_factor + def __init__(self, connection): self.connection = connection self.discr = connection.from_discr self.vis_discr = connection.to_discr @@ -139,8 +138,6 @@ class VTKVisualizer(object): def version(self): return "0.1" - # {{{ connectivity - @property def simplex_cell_types(self): import pyvisfile.vtk as vtk @@ -213,17 +210,9 @@ class VTKVisualizer(object): assert len(node_tuples) == grp.nunit_dofs return el_connectivity, vtk_cell_type - @memoize_method - def vis_nodes_numpy(self): - actx = self.vis_discr._setup_actx - return np.array([ - actx.to_numpy(flatten(thaw(actx, ary))) - for ary in self.vis_discr.nodes() - ]) - @property @memoize_method - def vtk_cells(self): + def cells(self): return np.hstack([ vgrp.vis_connectivity.reshape(-1) for vgrp in self.groups ]) @@ -261,94 +250,10 @@ class VTKVisualizer(object): return result - # }}} - - def write_vtk_file(self, file_name, names_and_fields, - compressor=None, real_only=False, overwrite=False): - from pyvisfile.vtk import ( - UnstructuredGrid, DataArray, - AppendedDataXMLGenerator, - VF_LIST_OF_COMPONENTS) - nodes = self.vis_nodes_numpy() - names_and_fields = [ - (name, resample_to_numpy(self.connection, fld)) - for name, fld in names_and_fields] +class VTKLagrangeConnectivity(VTKConnectivity): + """Connectivity for high-order Lagrange elements.""" - # {{{ create cell_types - - nsubelements = sum(vgrp.nsubelements for vgrp in self.groups) - cell_types = np.empty(nsubelements, dtype=np.uint8) - cell_types.fill(255) - - for vgrp in self.groups: - isubelements = np.s_[ - vgrp.subelement_nr_base: - vgrp.subelement_nr_base + vgrp.nsubelements] - cell_types[isubelements] = vgrp.vtk_cell_type - - assert (cell_types < 255).all() - - # }}} - - # {{{ shrink elements - - if abs(self.element_shrink_factor - 1.0) > 1.0e-14: - node_nr_base = 0 - for vgrp in self.vis_discr.groups: - nodes_view = ( - nodes[:, node_nr_base:node_nr_base + vgrp.ndofs] - .reshape(nodes.shape[0], vgrp.nelements, vgrp.nunit_dofs)) - - el_centers = np.mean(nodes_view, axis=-1) - nodes_view[:] = ( - (self.element_shrink_factor * nodes_view) - + (1-self.element_shrink_factor) - * el_centers[:, :, np.newaxis]) - - node_nr_base += vgrp.ndofs - - # }}} - - # {{{ create grid - - nodes = nodes.reshape(self.vis_discr.ambient_dim, -1) - points = DataArray("points", nodes, - vector_format=VF_LIST_OF_COMPONENTS) - - grid = UnstructuredGrid( - (nodes.shape[1], points), - cells=self.vtk_cells, - cell_types=cell_types) - - for name, field in separate_by_real_and_imag(names_and_fields, real_only): - grid.add_pointdata( - DataArray(name, field, vector_format=VF_LIST_OF_COMPONENTS) - ) - - # }}} - - # {{{ write - - import os - from meshmode import FileExistsError - if os.path.exists(file_name): - if overwrite: - os.remove(file_name) - else: - raise FileExistsError("output file '%s' already exists" % file_name) - - with open(file_name, "w") as outf: - generator = AppendedDataXMLGenerator( - compressor=compressor, - vtk_file_version=self.version) - - generator(grid).write(outf) - - # }}} - - -class VTKLagrangeVisualizer(VTKVisualizer): @property def version(self): # NOTE: version 2.2 has an updated ordering for the hexahedron @@ -382,7 +287,8 @@ class VTKLagrangeVisualizer(VTKVisualizer): vtk_lagrange_simplex_node_tuples, vtk_lagrange_simplex_node_tuples_to_permutation) - node_tuples = vtk_lagrange_simplex_node_tuples(grp.dim, grp.order) + node_tuples = vtk_lagrange_simplex_node_tuples( + grp.dim, grp.order, is_consistent=True) el_connectivity = np.array( vtk_lagrange_simplex_node_tuples_to_permutation(node_tuples), dtype=np.intp).reshape(1, 1, -1) @@ -398,7 +304,7 @@ class VTKLagrangeVisualizer(VTKVisualizer): @property @memoize_method - def vtk_cells(self): + def cells(self): connectivity = np.hstack([ grp.vis_connectivity.reshape(-1) for grp in self.groups @@ -436,18 +342,22 @@ class Visualizer(object): """ def __init__(self, connection, - element_shrink_factor=None, - use_high_order_vtk=False): + element_shrink_factor=None): self.connection = connection self.discr = connection.from_discr self.vis_discr = connection.to_discr - if use_high_order_vtk: - self.vtk = VTKLagrangeVisualizer(connection, - element_shrink_factor=element_shrink_factor) - else: - self.vtk = VTKVisualizer(connection, - element_shrink_factor=element_shrink_factor) + if element_shrink_factor is None: + element_shrink_factor = 1.0 + self.element_shrink_factor = element_shrink_factor + + @memoize_method + def _vis_nodes_numpy(self): + actx = self.vis_discr._setup_actx + return np.array([ + actx.to_numpy(flatten(thaw(actx, ary))) + for ary in self.vis_discr.nodes() + ]) # {{{ mayavi @@ -456,13 +366,11 @@ class Visualizer(object): do_show = kwargs.pop("do_show", True) - nodes = self.vtk.vis_nodes_numpy() + nodes = self._vis_nodes_numpy() field = resample_to_numpy(self.connection, field) assert nodes.shape[0] == self.vis_discr.ambient_dim - #mlab.points3d(nodes[0], nodes[1], 0*nodes[0]) - - vis_connectivity, = self._vis_connectivity() + vis_connectivity = self._vtk_connectivity.groups[0].vis_connectivity if self.vis_discr.dim == 1: nodes = list(nodes) @@ -502,14 +410,109 @@ class Visualizer(object): # {{{ vtk + @property + @memoize_method + def _vtk_connectivity(self): + return VTKConnectivity(self.connection) + + @property + @memoize_method + def _vtk_lagrange_connectivity(self): + return VTKLagrangeConnectivity(self.connection) + def write_vtk_file(self, file_name, names_and_fields, compressor=None, real_only=False, - overwrite=False): - self.vtk.write_vtk_file(file_name, names_and_fields, - compressor=compressor, - real_only=real_only, - overwrite=overwrite) + overwrite=False, + use_lagrange_elements=False): + from pyvisfile.vtk import ( + UnstructuredGrid, DataArray, + AppendedDataXMLGenerator, + VF_LIST_OF_COMPONENTS) + + if use_lagrange_elements: + connectivity = self._vtk_lagrange_connectivity + else: + connectivity = self._vtk_connectivity + + nodes = self._vis_nodes_numpy() + names_and_fields = [ + (name, resample_to_numpy(self.connection, fld)) + for name, fld in names_and_fields] + + # {{{ create cell_types + + nsubelements = sum(vgrp.nsubelements for vgrp in connectivity.groups) + cell_types = np.empty(nsubelements, dtype=np.uint8) + cell_types.fill(255) + + for vgrp in connectivity.groups: + isubelements = np.s_[ + vgrp.subelement_nr_base: + vgrp.subelement_nr_base + vgrp.nsubelements] + cell_types[isubelements] = vgrp.vtk_cell_type + + assert (cell_types < 255).all() + + # }}} + + # {{{ shrink elements + + if abs(self.element_shrink_factor - 1.0) > 1.0e-14: + node_nr_base = 0 + for vgrp in self.vis_discr.groups: + nodes_view = ( + nodes[:, node_nr_base:node_nr_base + vgrp.ndofs] + .reshape(nodes.shape[0], vgrp.nelements, vgrp.nunit_dofs)) + + el_centers = np.mean(nodes_view, axis=-1) + nodes_view[:] = ( + (self.element_shrink_factor * nodes_view) + + (1-self.element_shrink_factor) + * el_centers[:, :, np.newaxis]) + + node_nr_base += vgrp.ndofs + + # }}} + + # {{{ create grid + + nodes = nodes.reshape(self.vis_discr.ambient_dim, -1) + points = DataArray("points", nodes, + vector_format=VF_LIST_OF_COMPONENTS) + + grid = UnstructuredGrid( + (nodes.shape[1], points), + cells=connectivity.cells, + cell_types=cell_types) + + for name, field in separate_by_real_and_imag(names_and_fields, real_only): + grid.add_pointdata( + DataArray(name, field, vector_format=VF_LIST_OF_COMPONENTS) + ) + + # }}} + + # {{{ write + + import os + from meshmode import FileExistsError + if os.path.exists(file_name): + if overwrite: + os.remove(file_name) + else: + raise FileExistsError("output file '%s' already exists" % file_name) + + with open(file_name, "w") as outf: + generator = AppendedDataXMLGenerator( + compressor=compressor, + vtk_file_version=connectivity.version) + + generator(grid).write(outf) + + # }}} + + # }}} # {{{ matplotlib 3D @@ -524,12 +527,12 @@ class Visualizer(object): vmax = kwargs.pop("vmax", None) norm = kwargs.pop("norm", None) - nodes = self.vtk.vis_nodes_numpy() - field = resample_to_numpy(field) + nodes = self._vis_nodes_numpy() + field = resample_to_numpy(self.connection, field) assert nodes.shape[0] == self.vis_discr.ambient_dim - vis_connectivity, = self._vis_connectivity() + vis_connectivity, = self._vtk_connectivity.groups fig = plt.gcf() ax = fig.gca(projection="3d") @@ -580,18 +583,12 @@ class Visualizer(object): # }}} -def make_visualizer(actx, discr, vis_order, - element_shrink_factor=None, use_high_order_vtk=False): +def make_visualizer(actx, discr, vis_order, element_shrink_factor=None): from meshmode.discretization import Discretization - if use_high_order_vtk: - from meshmode.discretization.poly_element import ( - PolynomialEquidistantSimplexElementGroup as SimplexElementGroup, - EquidistantTensorProductElementGroup as TensorElementGroup) - else: - from meshmode.discretization.poly_element import ( - PolynomialWarpAndBlendElementGroup as SimplexElementGroup, - LegendreGaussLobattoTensorProductElementGroup as TensorElementGroup) + from meshmode.discretization.poly_element import ( + PolynomialEquidistantSimplexElementGroup as SimplexElementGroup, + EquidistantTensorProductElementGroup as TensorElementGroup) from meshmode.discretization.poly_element import OrderAndTypeBasedGroupFactory vis_discr = Discretization( @@ -607,8 +604,7 @@ def make_visualizer(actx, discr, vis_order, return Visualizer( make_same_mesh_connection(actx, vis_discr, discr), - element_shrink_factor=element_shrink_factor, - use_high_order_vtk=use_high_order_vtk) + element_shrink_factor=element_shrink_factor) # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index bfb6b539..cabbd593 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -84,6 +84,61 @@ def test_circle_mesh(visualize=False): # }}} +# {{{ test visualizer + +@pytest.mark.parametrize("dim", [1, 2, 3]) +def test_visualizers(ctx_factory, dim): + logging.basicConfig(level=logging.INFO) + + cl_ctx = ctx_factory() + queue = cl.CommandQueue(cl_ctx) + actx = PyOpenCLArrayContext(queue) + + nelements = 64 + target_order = 4 + + if dim == 1: + mesh = mgen.make_curve_mesh( + mgen.NArmedStarfish(5, 0.25), + np.linspace(0.0, 1.0, nelements + 1), + target_order) + elif dim == 2: + mesh = mgen.generate_torus(5.0, 1.0, order=target_order) + elif dim == 3: + mesh = mgen.generate_warped_rect_mesh(dim, target_order, 5) + else: + raise ValueError("unknown dimensionality") + + from meshmode.discretization import Discretization + discr = Discretization(actx, mesh, + InterpolatoryQuadratureSimplexGroupFactory(target_order)) + + from meshmode.discretization.visualization import make_visualizer + vis = make_visualizer(actx, discr, target_order) + + vis.write_vtk_file(f"visualizer_vtk_lagrange_{dim}.vtu", [], + use_lagrange_elements=True, overwrite=True) + vis.write_vtk_file(f"visualizer_vtk_linear_{dim}.vtu", [], + use_lagrange_elements=False, overwrite=True) + + if mesh.dim <= 2: + field = thaw(actx, discr.nodes()[0]) + + if mesh.dim == 2: + try: + vis.show_scalar_in_matplotlib_3d(field, do_show=False) + except ImportError: + logger.info("matplotlib not available") + + if mesh.dim <= 2: + try: + vis.show_scalar_in_mayavi(field, do_show=False) + except ImportError: + logger.info("mayavi not avaiable") + +# }}} + + # {{{ test boundary tags def test_boundary_tags(): -- GitLab From 19fa70778a0cd01af755d8799ca0ffc0e2d19a28 Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 17 Jul 2020 21:15:39 -0500 Subject: [PATCH 08/13] remove useless import as --- meshmode/discretization/visualization.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index 19f4f44c..b6319b3d 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -119,7 +119,7 @@ class _VisConnectivityGroup(Record): # }}} -# {{{ vtk visualizers +# {{{ vtk connectivity class VTKConnectivity: """Connectivity for standard linear VTK element types. @@ -341,8 +341,7 @@ class Visualizer(object): .. automethod:: write_vtk_file """ - def __init__(self, connection, - element_shrink_factor=None): + def __init__(self, connection, element_shrink_factor=None): self.connection = connection self.discr = connection.from_discr self.vis_discr = connection.to_discr @@ -478,8 +477,7 @@ class Visualizer(object): # {{{ create grid nodes = nodes.reshape(self.vis_discr.ambient_dim, -1) - points = DataArray("points", nodes, - vector_format=VF_LIST_OF_COMPONENTS) + points = DataArray("points", nodes, vector_format=VF_LIST_OF_COMPONENTS) grid = UnstructuredGrid( (nodes.shape[1], points), @@ -587,16 +585,17 @@ def make_visualizer(actx, discr, vis_order, element_shrink_factor=None): from meshmode.discretization import Discretization from meshmode.discretization.poly_element import ( - PolynomialEquidistantSimplexElementGroup as SimplexElementGroup, - EquidistantTensorProductElementGroup as TensorElementGroup) + PolynomialEquidistantSimplexElementGroup, + EquidistantTensorProductElementGroup, + OrderAndTypeBasedGroupFactory + ) - from meshmode.discretization.poly_element import OrderAndTypeBasedGroupFactory vis_discr = Discretization( actx, discr.mesh, OrderAndTypeBasedGroupFactory( vis_order, - simplex_group_class=SimplexElementGroup, - tensor_product_group_class=TensorElementGroup), + simplex_group_class=PolynomialEquidistantSimplexElementGroup, + tensor_product_group_class=EquidistantTensorProductElementGroup), real_dtype=discr.real_dtype) from meshmode.discretization.connection import \ -- GitLab From 89376f769ab6a36edffcc53a564e0236f2d695ed Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Fri, 17 Jul 2020 21:24:38 -0500 Subject: [PATCH 09/13] add pyvisfile to all the tests --- .gitlab-ci.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 3c9da790..b89185d7 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -2,7 +2,7 @@ Python 3 Nvidia K40: script: - export PY_EXE=python3 - export PYOPENCL_TEST=nvi:k40 - - export EXTRA_INSTALL="pybind11 cython numpy mako" + - export EXTRA_INSTALL="pybind11 cython numpy mako pyvisfile" # cython is here because pytential (for now, for TS) depends on it - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" @@ -20,7 +20,7 @@ Python 3 POCL: - export PY_EXE=python3 - export PYOPENCL_TEST=portable:pthread # cython is here because pytential (for now, for TS) depends on it - - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py" + - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py pyvisfile" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" tags: @@ -38,7 +38,7 @@ Python 3 POCL: - export PY_EXE=python3 - export PYOPENCL_TEST=portable:pthread # cython is here because pytential (for now, for TS) depends on it - - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py" + - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py pyvisfile" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" tags: -- GitLab From e128b67922e9827650e035950da282cfef12e9af Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 23 Jul 2020 11:10:26 -0500 Subject: [PATCH 10/13] add a write_high_order_vtk_file method to Visualizer --- meshmode/discretization/visualization.py | 65 +++++++++++++++++------- test/test_meshmode.py | 15 ++++-- 2 files changed, 57 insertions(+), 23 deletions(-) diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index b6319b3d..4a7c93f0 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -341,7 +341,8 @@ class Visualizer(object): .. automethod:: write_vtk_file """ - def __init__(self, connection, element_shrink_factor=None): + def __init__(self, connection, + element_shrink_factor=None, is_equidistant=False): self.connection = connection self.discr = connection.from_discr self.vis_discr = connection.to_discr @@ -349,6 +350,7 @@ class Visualizer(object): if element_shrink_factor is None: element_shrink_factor = 1.0 self.element_shrink_factor = element_shrink_factor + self.is_equidistant = is_equidistant @memoize_method def _vis_nodes_numpy(self): @@ -419,21 +421,33 @@ class Visualizer(object): def _vtk_lagrange_connectivity(self): return VTKLagrangeConnectivity(self.connection) + def write_high_order_vtk_file(self, file_name, names_and_fields, + compressor=None, real_only=False, overwrite=False): + if not self.is_equidistant: + raise RuntimeError("cannot visualize high-order Lagrange elements. " + "call 'make_visualizer' with 'force_equidistant=True'.") + + self._write_vtk_file(file_name, names_and_fields, + connectivity=self._vtk_lagrange_connectivity, + compressor=compressor, + real_only=real_only, + overwrite=overwrite) + def write_vtk_file(self, file_name, names_and_fields, - compressor=None, - real_only=False, - overwrite=False, - use_lagrange_elements=False): + compressor=None, real_only=False, overwrite=False): + self._write_vtk_file(file_name, names_and_fields, + connectivity=self._vtk_connectivity, + compressor=compressor, + real_only=real_only, + overwrite=overwrite) + + def _write_vtk_file(self, file_name, names_and_fields, connectivity, + compressor=None, real_only=False, overwrite=False): from pyvisfile.vtk import ( UnstructuredGrid, DataArray, AppendedDataXMLGenerator, VF_LIST_OF_COMPONENTS) - if use_lagrange_elements: - connectivity = self._vtk_lagrange_connectivity - else: - connectivity = self._vtk_connectivity - nodes = self._vis_nodes_numpy() names_and_fields = [ (name, resample_to_numpy(self.connection, fld)) @@ -581,21 +595,33 @@ class Visualizer(object): # }}} -def make_visualizer(actx, discr, vis_order, element_shrink_factor=None): +def make_visualizer(actx, discr, vis_order, + element_shrink_factor=None, force_equidistant=False): + """ + :arg vis_order: order of the visualization DOFs. + :arg element_shrink_factor: number in :math:`(0, 1]`. + :arg force_equidistant: if *True*, the visualization is done on + equidistant nodes. If plotting high-order Lagrange VTK elements, this + needs to be set to *True*. + """ from meshmode.discretization import Discretization - from meshmode.discretization.poly_element import ( - PolynomialEquidistantSimplexElementGroup, - EquidistantTensorProductElementGroup, - OrderAndTypeBasedGroupFactory - ) + if force_equidistant: + from meshmode.discretization.poly_element import ( + PolynomialEquidistantSimplexElementGroup as SimplexElementGroup, + EquidistantTensorProductElementGroup as TensorElementGroup) + else: + from meshmode.discretization.poly_element import ( + PolynomialWarpAndBlendElementGroup as SimplexElementGroup, + LegendreGaussLobattoTensorProductElementGroup as TensorElementGroup) + from meshmode.discretization.poly_element import OrderAndTypeBasedGroupFactory vis_discr = Discretization( actx, discr.mesh, OrderAndTypeBasedGroupFactory( vis_order, - simplex_group_class=PolynomialEquidistantSimplexElementGroup, - tensor_product_group_class=EquidistantTensorProductElementGroup), + simplex_group_class=SimplexElementGroup, + tensor_product_group_class=TensorElementGroup), real_dtype=discr.real_dtype) from meshmode.discretization.connection import \ @@ -603,7 +629,8 @@ def make_visualizer(actx, discr, vis_order, element_shrink_factor=None): return Visualizer( make_same_mesh_connection(actx, vis_discr, discr), - element_shrink_factor=element_shrink_factor) + element_shrink_factor=element_shrink_factor, + is_equidistant=force_equidistant) # }}} diff --git a/test/test_meshmode.py b/test/test_meshmode.py index cabbd593..48946bce 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -116,10 +116,12 @@ def test_visualizers(ctx_factory, dim): from meshmode.discretization.visualization import make_visualizer vis = make_visualizer(actx, discr, target_order) - vis.write_vtk_file(f"visualizer_vtk_lagrange_{dim}.vtu", [], - use_lagrange_elements=True, overwrite=True) - vis.write_vtk_file(f"visualizer_vtk_linear_{dim}.vtu", [], - use_lagrange_elements=False, overwrite=True) + vis.write_vtk_file(f"visualizer_vtk_linear_{dim}.vtu", + [], overwrite=True) + + with pytest.raises(RuntimeError): + vis.write_high_order_vtk_file(f"visualizer_vtk_lagrange_{dim}.vtu", + [], overwrite=True) if mesh.dim <= 2: field = thaw(actx, discr.nodes()[0]) @@ -136,6 +138,11 @@ def test_visualizers(ctx_factory, dim): except ImportError: logger.info("mayavi not avaiable") + vis = make_visualizer(actx, discr, target_order, + force_equidistant=True) + vis.write_high_order_vtk_file(f"visualizer_vtk_lagrange_{dim}.vtu", + [], overwrite=True) + # }}} -- GitLab From 7b36d7dddbb8b077acc684637c68fdc57ed83d1e Mon Sep 17 00:00:00 2001 From: Alexandru Fikl Date: Thu, 23 Jul 2020 11:16:49 -0500 Subject: [PATCH 11/13] add a docstring to write_high_order_vtk_file --- meshmode/discretization/visualization.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index 4a7c93f0..1ef1d121 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -339,6 +339,7 @@ class Visualizer(object): .. automethod:: show_scalar_in_mayavi .. automethod:: show_scalar_in_matplotlib_3d .. automethod:: write_vtk_file + .. automethod:: write_high_order_vtk_file """ def __init__(self, connection, @@ -423,6 +424,10 @@ class Visualizer(object): def write_high_order_vtk_file(self, file_name, names_and_fields, compressor=None, real_only=False, overwrite=False): + """Writes arbitrary order Lagrange VTK elements. These elements are + described in `this blog post `_ + and are available in VTK 8.1 and newer. + """ # noqa if not self.is_equidistant: raise RuntimeError("cannot visualize high-order Lagrange elements. " "call 'make_visualizer' with 'force_equidistant=True'.") -- GitLab From 86c71fff0406587db74c4d53d976502f12741f32 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 28 Jul 2020 17:31:41 -0500 Subject: [PATCH 12/13] Move pyvisfile from CI config to requirements.txt --- .gitlab-ci.yml | 6 +++--- requirements.txt | 1 + 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index b89185d7..3c9da790 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -2,7 +2,7 @@ Python 3 Nvidia K40: script: - export PY_EXE=python3 - export PYOPENCL_TEST=nvi:k40 - - export EXTRA_INSTALL="pybind11 cython numpy mako pyvisfile" + - export EXTRA_INSTALL="pybind11 cython numpy mako" # cython is here because pytential (for now, for TS) depends on it - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" @@ -20,7 +20,7 @@ Python 3 POCL: - export PY_EXE=python3 - export PYOPENCL_TEST=portable:pthread # cython is here because pytential (for now, for TS) depends on it - - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py pyvisfile" + - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" tags: @@ -38,7 +38,7 @@ Python 3 POCL: - export PY_EXE=python3 - export PYOPENCL_TEST=portable:pthread # cython is here because pytential (for now, for TS) depends on it - - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py pyvisfile" + - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py" - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh - ". ./build-and-test-py-project.sh" tags: diff --git a/requirements.txt b/requirements.txt index 574e5094..939b0e4b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,6 +2,7 @@ numpy recursivenodes git+https://gitlab.tiker.net/inducer/pytools.git#egg=pytools git+https://gitlab.tiker.net/inducer/gmsh_interop.git#egg=gmsh_interop +git+https://gitlab.tiker.net/inducer/pyvisfile.git#egg=pyvisfile git+https://gitlab.tiker.net/inducer/modepy.git#egg=modepy git+https://gitlab.tiker.net/inducer/pyopencl.git#egg=pyopencl git+https://gitlab.tiker.net/inducer/islpy.git#egg=islpy -- GitLab From 6e266de37418152955aabf3a767eab374c3164c9 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 28 Jul 2020 17:32:16 -0500 Subject: [PATCH 13/13] Point requirements.txt URLs at Github instead of Gitlab --- requirements.txt | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/requirements.txt b/requirements.txt index 939b0e4b..6d3fc187 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,18 +1,18 @@ numpy recursivenodes -git+https://gitlab.tiker.net/inducer/pytools.git#egg=pytools -git+https://gitlab.tiker.net/inducer/gmsh_interop.git#egg=gmsh_interop -git+https://gitlab.tiker.net/inducer/pyvisfile.git#egg=pyvisfile -git+https://gitlab.tiker.net/inducer/modepy.git#egg=modepy -git+https://gitlab.tiker.net/inducer/pyopencl.git#egg=pyopencl -git+https://gitlab.tiker.net/inducer/islpy.git#egg=islpy +git+https://github.com/inducer/pytools.git#egg=pytools +git+https://github.com/inducer/gmsh_interop.git#egg=gmsh_interop +git+https://github.com/inducer/pyvisfile.git#egg=pyvisfile +git+https://github.com/inducer/modepy.git#egg=modepy +git+https://github.com/inducer/pyopencl.git#egg=pyopencl +git+https://github.com/inducer/islpy.git#egg=islpy # required by pytential, which is in turn needed for some tests -git+https://gitlab.tiker.net/inducer/pymbolic.git#egg=pymbolic +git+https://github.com/inducer/pymbolic.git#egg=pymbolic # also depends on pymbolic, so should come after it -git+https://gitlab.tiker.net/inducer/loopy.git#egg=loo.py +git+https://github.com/inducer/loopy.git#egg=loo.py # more pytential dependencies git+https://github.com/inducer/boxtree.git#egg=boxtree @@ -20,4 +20,4 @@ git+https://github.com/inducer/sumpy.git#egg=sumpy git+https://github.com/inducer/pytential.git#pytential # requires pymetis for tests for partition_mesh -git+https://gitlab.tiker.net/inducer/pymetis.git#egg=pymetis +git+https://github.com/inducer/pymetis.git#egg=pymetis -- GitLab