diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py
index 7129c1dcf4b7268d4bc62dcccaab8cae6268dc95..d4d07ac4903308ed9ea371fc19e1b4ccdfa5073b 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
^^^^^^^^^^^^^^^
@@ -342,16 +343,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..5e20d2216aea16211fcdc7aba6b109d4bed24c05 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,107 @@ 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):
+ assert self.is_equidistant
+ 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 <https://blog.kitware.com/modeling-arbitrary-order-lagrange-finite-elements-in-the-visualization-toolkit/>`_
+ and are available in VTK 8.1 and newer.
+ """ # noqa
+ if not self.is_equidistant:
+ raise RuntimeError("Cannot visualize high-order Lagrange elements "
+ "using a non-equidistant visualizer. "
+ "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 +523,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 +546,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 +602,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)