diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index e317e31130237ad2adf5173af4a7c0d23d0f2b9d..55cce2eab1c9b994225b5ea8caa91dfa21912523 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -43,6 +43,7 @@ Group types .. autoclass:: PolynomialRecursiveNodesElementGroup .. autoclass:: PolynomialEquidistantSimplexElementGroup .. autoclass:: LegendreGaussLobattoTensorProductElementGroup +.. autoclass:: EquidistantTensorProductElementGroup Group factories ^^^^^^^^^^^^^^^ @@ -310,16 +311,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 24fe92f675293ef6ab048553778ac50a394d4173..1ef1d1216ed45e512fab5b4383d36676ea3d9771 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -35,13 +35,15 @@ __doc__ = """ """ -# {{{ visualizer +# {{{ helpers -def separate_by_real_and_imag(data, real_only): - # This function is called on numpy data that has already been - # merged into a single vector. +def separate_by_real_and_imag(names_and_fields, real_only): + """ + :arg names_and_fields: input data array must be already flattened into a + single :mod:`numpy` array using :func:`resample_to_numpy`. + """ - for name, field in data: + for name, field in names_and_fields: if isinstance(field, np.ndarray) and field.dtype.char == "O": assert len(field.shape) == 1 from pytools.obj_array import ( @@ -53,25 +55,45 @@ def separate_by_real_and_imag(data, real_only): yield (name, obj_array_vectorize(obj_array_real_copy, field)) else: - yield (name+"_r", + yield (f"{name}_r", obj_array_vectorize(obj_array_real_copy, field)) - yield (name+"_i", + yield (f"{name}_i", obj_array_vectorize(obj_array_imag_copy, field)) else: yield (name, field) else: if field.dtype.kind == "c": - yield (name+"_r", field.real.copy()) - yield (name+"_i", field.imag.copy()) + if real_only: + yield (name, field.real.copy()) + else: + yield (f"{name}_r", field.real.copy()) + yield (f"{name}_i", field.imag.copy()) else: yield (name, field) +def resample_to_numpy(conn, vec): + if (isinstance(vec, np.ndarray) + and vec.dtype.char == "O" + and not isinstance(vec, DOFArray)): + from pytools.obj_array import obj_array_vectorize + return obj_array_vectorize(lambda x: resample_to_numpy(conn, x), vec) + + from numbers import Number + if isinstance(vec, Number): + nnodes = sum(grp.ndofs for grp in conn.to_discr.groups) + return np.ones(nnodes) * vec + else: + resampled = conn(vec) + actx = resampled.array_context + return actx.to_numpy(flatten(resampled)) + + 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 @@ -94,130 +116,128 @@ class _VisConnectivityGroup(Record): def primitive_element_size(self): return self.vis_connectivity.shape[2] +# }}} -class Visualizer(object): - """ - .. automethod:: show_scalar_in_mayavi - .. automethod:: show_scalar_in_matplotlib_3d - .. automethod:: write_vtk_file + +# {{{ vtk connectivity + +class VTKConnectivity: + """Connectivity for standard linear VTK element types. + + .. attribute:: version + .. attribute:: cells + .. attribute:: groups """ - def __init__(self, connection, element_shrink_factor=None): + def __init__(self, connection): 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 + @property + def simplex_cell_types(self): + import pyvisfile.vtk as vtk + return { + 1: vtk.VTK_LINE, + 2: vtk.VTK_TRIANGLE, + 3: vtk.VTK_TETRA, + } + + @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 simplex_submesh + el_connectivity = np.array( + simplex_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] - def _resample_to_numpy(self, vec): - if (isinstance(vec, np.ndarray) - and vec.dtype.char == "O" - and not isinstance(vec, DOFArray)): - from pytools.obj_array import obj_array_vectorize - return obj_array_vectorize(self._resample_to_numpy, vec) - - from numbers import Number - if isinstance(vec, Number): - nnodes = sum(grp.ndofs for grp in self.connection.to_discr.groups) - return np.ones(nnodes) * vec else: - resampled = self.connection(vec) - actx = resampled.array_context - return actx.to_numpy(flatten(resampled)) + raise NotImplementedError("visualization for element groups " + "of type '%s'" % type(grp.mesh_el_group).__name__) + + assert len(node_tuples) == grp.nunit_dofs + return el_connectivity, vtk_cell_type + @property @memoize_method - def _vis_nodes(self): - actx = self.vis_discr._setup_actx - return np.array([ - actx.to_numpy(flatten(thaw(actx, ary))) - for ary in self.vis_discr.nodes() + def cells(self): + return np.hstack([ + vgrp.vis_connectivity.reshape(-1) for vgrp in self.groups ]) - # {{{ vis sub-element connectivity - + @property @memoize_method - def _vis_connectivity(self): + 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 node_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 simplex_submesh - el_connectivity = np.array( - simplex_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] - - else: - raise NotImplementedError("visualization for element groups " - "of type '%s'" % type(group.mesh_el_group).__name__) + for grp in self.vis_discr.groups: + el_connectivity, vtk_cell_type = \ + self.connectivity_for_element_group(grp) - assert len(node_tuples) == group.nunit_dofs - vis_connectivity = ( - node_nr_base + np.arange( - 0, group.nelements*group.nunit_dofs, group.nunit_dofs - )[:, np.newaxis, np.newaxis] - + el_connectivity).astype(np.intp) + offsets = node_nr_base + np.arange( + 0, + grp.nelements * grp.nunit_dofs, + grp.nunit_dofs).reshape(-1, 1, 1) + vis_connectivity = (offsets + el_connectivity).astype(np.intp) vgrp = _VisConnectivityGroup( vis_connectivity=vis_connectivity, @@ -226,11 +246,120 @@ class Visualizer(object): result.append(vgrp) subel_nr_base += vgrp.nsubelements - node_nr_base += group.ndofs + node_nr_base += grp.ndofs return result - # }}} + +class VTKLagrangeConnectivity(VTKConnectivity): + """Connectivity for high-order Lagrange elements.""" + + @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.0" + + @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, is_consistent=True) + 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_dofs + 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 + ]) + + grp_offsets = np.cumsum([0] + [ + grp.ndofs for grp in self.vis_discr.groups + ]) + + offsets = np.hstack([ + grp_offset + np.arange( + grp.nunit_dofs, + grp.nelements * grp.nunit_dofs + 1, + grp.nunit_dofs) + for grp_offset, grp in zip(grp_offsets, 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 + .. automethod:: write_high_order_vtk_file + """ + + 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 + + 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): + actx = self.vis_discr._setup_actx + return np.array([ + actx.to_numpy(flatten(thaw(actx, ary))) + for ary in self.vis_discr.nodes() + ]) # {{{ mayavi @@ -239,13 +368,11 @@ class Visualizer(object): do_show = kwargs.pop("do_show", True) - nodes = self._vis_nodes() - field = self._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 - #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) @@ -285,67 +412,105 @@ class Visualizer(object): # {{{ vtk - def write_vtk_file(self, file_name, names_and_fields, - compressor=None, - real_only=False, - overwrite=False): + @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_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'.") + + 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): + 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) - nodes = self._vis_nodes() + nodes = self._vis_nodes_numpy() names_and_fields = [ - (name, self._resample_to_numpy(fld)) + (name, resample_to_numpy(self.connection, 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) + nsubelements = sum(vgrp.nsubelements for vgrp in connectivity.groups) cell_types = np.empty(nsubelements, dtype=np.uint8) cell_types.fill(255) - for vgrp in vc_groups: - cell_types[ + + for vgrp in connectivity.groups: + isubelements = np.s_[ vgrp.subelement_nr_base: - vgrp.subelement_nr_base + vgrp.nsubelements] = \ - vgrp.vtk_cell_type + vgrp.subelement_nr_base + vgrp.nsubelements] + cell_types[isubelements] = vgrp.vtk_cell_type + assert (cell_types < 255).all() # }}} - if self.element_shrink_factor != 1: + # {{{ 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 = vgrp.view(nodes) + 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]) - if len(self.vis_discr.groups) != 1: - raise NotImplementedError("visualization with multiple " - "element groups") + 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( - (sum(grp.ndofs for grp in self.vis_discr.groups), - 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]), + (nodes.shape[1], points), + cells=connectivity.cells, 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)) + grid.add_pointdata( + DataArray(name, field, vector_format=VF_LIST_OF_COMPONENTS) + ) + + # }}} + + # {{{ write import os from meshmode import FileExistsError @@ -356,10 +521,16 @@ class Visualizer(object): raise FileExistsError("output file '%s' already exists" % file_name) with open(file_name, "w") as outf: - AppendedDataXMLGenerator(compressor)(grid).write(outf) + generator = AppendedDataXMLGenerator( + compressor=compressor, + vtk_file_version=connectivity.version) + + generator(grid).write(outf) # }}} + # }}} + # {{{ matplotlib 3D def show_scalar_in_matplotlib_3d(self, field, **kwargs): @@ -373,12 +544,12 @@ class Visualizer(object): vmax = kwargs.pop("vmax", None) norm = kwargs.pop("norm", None) - nodes = self._vis_nodes() - field = self._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") @@ -429,26 +600,42 @@ 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 ( - PolynomialWarpAndBlendElementGroup, - LegendreGaussLobattoTensorProductElementGroup, - 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=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(actx, vis_discr, discr), - element_shrink_factor=element_shrink_factor) + element_shrink_factor=element_shrink_factor, + is_equidistant=force_equidistant) # }}} diff --git a/requirements.txt b/requirements.txt index 574e50949645b761b6f3579c87fbb4e1dbfecc5e..6d3fc18752d7f8d2e4a3eb04a80f85b61b050e32 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,17 +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/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 @@ -19,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 diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 3417d0aa5128f7df3eb22f5c912a675dfabcb79e..48946bceb665c95d41c59a3f035147197d14adca 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,84 @@ 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) + +# }}} + + +# {{{ 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_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]) + + 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") + + vis = make_visualizer(actx, discr, target_order, + force_equidistant=True) + vis.write_high_order_vtk_file(f"visualizer_vtk_lagrange_{dim}.vtu", + [], overwrite=True) # }}} @@ -558,7 +624,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) @@ -573,7 +639,7 @@ def test_3d_orientation(ctx_factory, what, mesh_gen_func, visualize=False): from meshmode.discretization import Discretization discr = Discretization(actx, mesh, - PolynomialWarpAndBlendGroupFactory(1)) + PolynomialWarpAndBlendGroupFactory(3)) from pytential import bind, sym @@ -604,9 +670,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(actx, discr, 1) + vis = make_visualizer(actx, discr, 3) - vis.write_vtk_file("normals.vtu", [ + vis.write_vtk_file("orientation_3d_%s_normals.vtu" % what, [ ("normals", normals), ]) @@ -647,7 +713,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 @@ -658,12 +724,12 @@ 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.write_vtk_file("merged.vtu", []) + vis = make_visualizer(actx, discr, 3, element_shrink_factor=0.8) + vis.write_vtk_file("merge_and_map.vtu", []) # }}} @@ -728,20 +794,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(actx, bdry_discr, 4) - - # }}} - from pytential import bind, sym bdry_normals = bind(bdry_discr, sym.normal(dim))(actx).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(actx, bdry_discr, 4) + + bdry_vis.write_vtk_file("sanity_single_element_boundary.vtu", [ + ("normals", bdry_normals) ]) normal_outward_check = bind(bdry_discr, @@ -858,14 +919,6 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, visualize=False): # }}} - # {{{ visualizers - - from meshmode.discretization.visualization import make_visualizer - vol_vis = make_visualizer(actx, vol_discr, 20) - bdry_vis = make_visualizer(actx, bdry_discr, 20) - - # }}} - from math import gamma true_surf = 2*np.pi**(dim/2)/gamma(dim/2) true_vol = true_surf/dim @@ -894,14 +947,22 @@ def test_sanity_balls(ctx_factory, src_file, dim, mesh_order, visualize=False): 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(actx, vol_discr, 7) + bdry_vis = make_visualizer(actx, 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)) (actx)), ]) - 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 @@ -931,11 +992,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 @@ -952,13 +1013,14 @@ 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.write_vtk_file("box.vtu", []) + vis = make_visualizer(actx, discr, 7) + vis.write_vtk_file("box_mesh.vtu", []) # }}} @@ -995,7 +1057,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) @@ -1013,7 +1075,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) @@ -1159,7 +1221,7 @@ def test_vtk_overwrite(ctx_factory): import os from meshmode import FileExistsError - filename = "test_vtk_overwrite.vtu" + filename = "vtk_overwrite_temp.vtu" if os.path.exists(filename): os.remove(filename) @@ -1410,13 +1472,12 @@ 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) - vis.write_vtk_file("test_mesh_multiple_groups.vtu", [ + vis.write_vtk_file("mesh_multiple_groups.vtu", [ ("group_id", group_id) ], overwrite=True)