diff --git a/.editorconfig b/.editorconfig
new file mode 100644
index 0000000000000000000000000000000000000000..dcbc21d86f9e4b17ea7e8803d538c4c0f0b6276a
--- /dev/null
+++ b/.editorconfig
@@ -0,0 +1,32 @@
+# https://editorconfig.org/
+# https://github.com/editorconfig/editorconfig-vim
+# https://github.com/editorconfig/editorconfig-emacs
+
+root = true
+
+[*]
+indent_style = space
+end_of_line = lf
+charset = utf-8
+trim_trailing_whitespace = true
+insert_final_newline = true
+
+[*.py]
+indent_size = 4
+
+[*.rst]
+indent_size = 4
+
+[*.cpp]
+indent_size = 2
+
+[*.hpp]
+indent_size = 2
+
+# There may be one in doc/
+[Makefile]
+indent_style = tab
+
+# https://github.com/microsoft/vscode/issues/1679
+[*.md]
+trim_trailing_whitespace = false
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 01e2cbfc5ee45291dbc4054c48ef98aec487e150..1e93235eb01aadc42a04c867974cc98ebf3ad143 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -24,19 +24,6 @@ jobs:
curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/prepare-and-run-flake8.sh
. ./prepare-and-run-flake8.sh "$(basename $GITHUB_REPOSITORY)" ./test
- pytest2:
- name: Pytest Conda Py2
- runs-on: ubuntu-latest
- steps:
- - uses: actions/checkout@v2
- - name: "Main Script"
- run: |
- sed 's/python=3/python=2.7/' .test-conda-env-py3.yml > .test-conda-env-py2.yml
- cat .test-conda-env-py2.yml
- CONDA_ENVIRONMENT=.test-conda-env-py2.yml
- curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project-within-miniconda.sh
- . ./build-and-test-py-project-within-miniconda.sh
-
pytest3:
name: Pytest Conda Py3
runs-on: ubuntu-latest
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 025048aca2f6cc88315050131c7b24801f1841d6..195145687eee9bb3495176547933715e801426a6 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,21 +1,3 @@
-Python 2.7 POCL:
- script:
- - export PY_EXE=python2.7
- - export PYOPENCL_TEST=portable
- - export EXTRA_INSTALL="pybind11 cython numpy mako mpi4py"
- # 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"
- tags:
- - python2.7
- - pocl
- - mpi
- except:
- - tags
- artifacts:
- reports:
- junit: test/pytest.xml
-
Python 3 Nvidia K40:
script:
- export PY_EXE=python3
@@ -36,7 +18,7 @@ Python 3 Nvidia K40:
Python 3 POCL:
script:
- export PY_EXE=python3
- - export PYOPENCL_TEST=portable
+ - 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"
- curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh
@@ -54,7 +36,7 @@ Python 3 POCL:
Python 3 POCL:
script:
- export PY_EXE=python3
- - export PYOPENCL_TEST=portable
+ - 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"
- curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-and-test-py-project.sh
diff --git a/doc/misc.rst b/doc/misc.rst
index 488262acbc29882a93c25d4ba137d77ac37e8660..b3fa9ffae0c88eef24b42f1109e463b0ab2021ea 100644
--- a/doc/misc.rst
+++ b/doc/misc.rst
@@ -3,28 +3,44 @@
Installation
============
-This command should install :mod:`meshmode`::
+This set of instructions is intended for 64-bit Linux and MacOS computers.
- pip install meshmode
+#. Make sure your system has the basics to build software.
-(Note the extra "."!)
+ On Debian derivatives (Ubuntu and many more),
+ installing ``build-essential`` should do the trick.
-You may need to run this with :command:`sudo`.
-If you don't already have `pip <https://pypi.python.org/pypi/pip>`_,
-run this beforehand::
+ Everywhere else, just making sure you have the ``g++`` package should be
+ enough.
- curl -O https://raw.github.com/pypa/pip/master/contrib/get-pip.py
- python get-pip.py
+#. Installing `miniforge for Python 3 on your respective system <https://github.com/conda-forge/miniforge>`_.
-For a more manual installation, `download the source
-<http://pypi.python.org/pypi/meshmode>`_, unpack it, and say::
+#. ``export CONDA=/WHERE/YOU/INSTALLED/miniforge3``
- python setup.py install
+ If you accepted the default location, this should work:
-You may also clone its git repository::
+ ``export CONDA=$HOME/miniforge3``
- git clone --recursive git://github.com/inducer/meshmode
- git clone --recursive http://git.tiker.net/trees/meshmode.git
+#. ``$CONDA/bin/conda create -n dgfem``
+
+#. ``source $CONDA/bin/activate dgfem``
+
+#. ``conda config --add channels conda-forge``
+
+#. ``conda install git pip pocl islpy pyopencl``
+
+#. Type the following command::
+
+ hash -r; for i in pymbolic cgen genpy modepy pyvisfile loopy meshmode; do python -m pip install git+https://github.com/inducer/$i.git; done
+
+Next time you want to use :mod:`meshmode`, just run the following command::
+
+ source /WHERE/YOU/INSTALLED/miniforge3/bin/activate dgfem
+
+You may also like to add this to a startup file (like :file:`$HOME/.bashrc`) or create an alias for it.
+
+After this, you should be able to run the `tests <https://github.com/inducer/meshmode/tree/master/test>`_
+or `examples <https://github.com/inducer/meshmode/tree/master/examples>`_.
User-visible Changes
====================
diff --git a/meshmode/discretization/connection/opposite_face.py b/meshmode/discretization/connection/opposite_face.py
index 97570b338f9938165004a46c59fb3ec1570e0cbe..e084525ba9a6eb5001548f432c25a0a8f861173e 100644
--- a/meshmode/discretization/connection/opposite_face.py
+++ b/meshmode/discretization/connection/opposite_face.py
@@ -342,22 +342,42 @@ def make_opposite_face_connection(volume_to_bdry_conn):
# {{{ index wrangling
- # Assert that the adjacency group and the restriction
- # interpolation batch and the adjacency group have the same
- # element ordering.
+ # The elements in the adjacency group will be a subset of
+ # the elements in the restriction interpolation batch:
+ # Imagine an inter-group boundary. The volume-to-boundary
+ # connection will include all faces as targets, whereas
+ # there will be separate adjacency groups for intra- and
+ # inter-group connections.
adj_tgt_flags = adj.element_faces == i_face_tgt
+ adj_els = adj.elements[adj_tgt_flags]
+ if adj_els.size == 0:
+ # NOTE: this case can happen for inter-group boundaries
+ # when all elements are adjacent on the same face
+ # index, so all other ones will be empty
+ continue
+
+ vbc_els = vbc_tgt_grp_face_batch.from_element_indices.get(queue)
+
+ if len(adj_els) == len(vbc_els):
+ # Same length: assert (below) that the two use the same
+ # ordering.
+ vbc_used_els = slice(None)
+
+ else:
+ # Genuine subset: figure out an index mapping.
+ vbc_els_sort_idx = np.argsort(vbc_els)
+ vbc_used_els = vbc_els_sort_idx[np.searchsorted(
+ vbc_els, adj_els, sorter=vbc_els_sort_idx
+ )]
- assert (np.array_equal(
- adj.elements[adj_tgt_flags],
- vbc_tgt_grp_face_batch.from_element_indices
- .get(queue=queue)))
+ assert np.array_equal(vbc_els[vbc_used_els], adj_els)
# find to_element_indices
tgt_bdry_element_indices = (
vbc_tgt_grp_face_batch.to_element_indices
- .get(queue=queue))
+ .get(queue=queue)[vbc_used_els])
# find from_element_indices
diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py
index 9df70d387b5d40d435982024c1dd0b78429915d8..c2ace66b3e8f1d97f3b98f0d0538d6deded87c75 100644
--- a/meshmode/mesh/__init__.py
+++ b/meshmode/mesh/__init__.py
@@ -1382,7 +1382,7 @@ def is_boundary_tag_empty(mesh, boundary_tag):
# }}}
-# {{{
+# {{{ is_affine_simplex_group
def is_affine_simplex_group(group, abs_tol=None):
if abs_tol is None:
diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py
index 8e5288bb8e05b7d27a6ed68d0cb364a546ecf11c..4fa520bc4d7a62a4dfb27f55025ff9d906b52a0f 100644
--- a/meshmode/mesh/generation.py
+++ b/meshmode/mesh/generation.py
@@ -446,7 +446,7 @@ def generate_icosphere(r, order, uniform_refinement_rounds=0):
mesh = generate_icosahedron(r, order)
if uniform_refinement_rounds:
- # These come out conformal, so we're OK to use the faster refiner.
+ # These come out conforming, so we're OK to use the faster refiner.
from meshmode.mesh.refinement import RefinerWithoutAdjacency
refiner = RefinerWithoutAdjacency(mesh)
for i in range(uniform_refinement_rounds):
@@ -673,7 +673,7 @@ def generate_urchin(order, m, n, est_rel_interp_tolerance, min_rad=0.2):
# {{{ generate_box_mesh
def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
- group_factory=None):
+ group_factory=None, boundary_tag_to_face=None):
"""Create a semi-structured mesh.
:param axis_coords: a tuple with a number of entries corresponding
@@ -681,12 +681,27 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
specifying the coordinates to be used along that axis.
:param group_factory: One of :class:`meshmode.mesh.SimplexElementGroup`
or :class:`meshmode.mesh.TensorProductElementGroup`.
+ :param boundary_tag_to_face: an optional dictionary for tagging boundaries.
+ The keys correspond to custom boundary tags, with the values giving
+ a list of the faces on which they should be applied in terms of coordinate
+ directions (``+x``, ``-x``, ``+y``, ``-y``, ``+z``, ``-z``, ``+w``, ``-w``).
+
+ For example::
+
+ boundary_tag_to_face={"bdry_1": ["+x", "+y"], "bdry_2": ["-x"]}
.. versionchanged:: 2017.1
*group_factory* parameter added.
+
+ .. versionchanged:: 2020.1
+
+ *boundary_tag_to_face* parameter added.
"""
+ if boundary_tag_to_face is None:
+ boundary_tag_to_face = {}
+
for iaxis, axc in enumerate(axis_coords):
if len(axc) < 2:
raise ValueError("need at least two points along axis %d"
@@ -780,8 +795,7 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
el_vertices.append((a011, a111, a101, a110))
else:
- raise NotImplementedError("box meshes of dimension %d"
- % dim)
+ raise NotImplementedError("box meshes of dimension %d" % dim)
el_vertices = np.array(el_vertices, dtype=np.int32)
@@ -789,22 +803,85 @@ def generate_box_mesh(axis_coords, order=1, coord_dtype=np.float64,
vertices.reshape(dim, -1), el_vertices, order,
group_factory=group_factory)
+ # {{{ compute facial adjacency for mesh if there is tag information
+
+ facial_adjacency_groups = None
+ face_vertex_indices_to_tags = {}
+ boundary_tags = list(boundary_tag_to_face.keys())
+ axes = ["x", "y", "z", "w"]
+
+ if boundary_tags:
+ from pytools import indices_in_shape
+ vert_index_to_tuple = {
+ vertex_indices[itup]: itup
+ for itup in indices_in_shape(shape)}
+
+ for tag_idx, tag in enumerate(boundary_tags):
+ # Need to map the correct face vertices to the boundary tags
+ for face in boundary_tag_to_face[tag]:
+ if len(face) != 2:
+ raise ValueError("face identifier '%s' does not "
+ "consist of exactly two characters" % face)
+
+ side, axis = face
+ try:
+ axis = axes.index(axis)
+ except ValueError:
+ raise ValueError("unrecognized axis in face identifier '%s'" % face)
+ if axis >= dim:
+ raise ValueError("axis in face identifier '%s' does not exist in %dD"
+ % (face, dim))
+
+ if side == "-":
+ vert_crit = 0
+ elif side == "+":
+ vert_crit = shape[axis] - 1
+ else:
+ raise ValueError("first character of face identifier '%s' is not"
+ "'+' or '-'" % face)
+
+ for ielem in range(0, grp.nelements):
+ for ref_fvi in grp.face_vertex_indices():
+ fvi = grp.vertex_indices[ielem, ref_fvi]
+ fvi_tuples = [vert_index_to_tuple[i] for i in fvi]
+
+ if all(fvi_tuple[axis] == vert_crit for fvi_tuple in fvi_tuples):
+ key = frozenset(fvi)
+ face_vertex_indices_to_tags.setdefault(key, []).append(tag)
+
+ if boundary_tags:
+ from meshmode.mesh import (
+ _compute_facial_adjacency_from_vertices, BTAG_ALL, BTAG_REALLY_ALL)
+ boundary_tags.extend([BTAG_ALL, BTAG_REALLY_ALL])
+ facial_adjacency_groups = _compute_facial_adjacency_from_vertices(
+ [grp], boundary_tags, np.int32, np.int8,
+ face_vertex_indices_to_tags)
+ else:
+ facial_adjacency_groups = None
+
+ # }}}
+
from meshmode.mesh import Mesh
return Mesh(vertices, [grp],
- is_conforming=True)
+ facial_adjacency_groups=facial_adjacency_groups,
+ is_conforming=True, boundary_tags=boundary_tags)
# }}}
# {{{ generate_regular_rect_mesh
-def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
+def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1,
+ boundary_tag_to_face=None,
+ group_factory=None):
"""Create a semi-structured rectangular mesh.
:param a: the lower left hand point of the rectangle
:param b: the upper right hand point of the rectangle
:param n: a tuple of integers indicating the total number of points
on [a,b].
+ :param boundary_tag_to_face: an optional dictionary for tagging boundaries.
+ See :func:`generate_box_mesh`.
"""
if min(n) < 2:
raise ValueError("need at least two points in each direction")
@@ -812,7 +889,9 @@ def generate_regular_rect_mesh(a=(0, 0), b=(1, 1), n=(5, 5), order=1):
axis_coords = [np.linspace(a_i, b_i, n_i)
for a_i, b_i, n_i in zip(a, b, n)]
- return generate_box_mesh(axis_coords, order=order)
+ return generate_box_mesh(axis_coords, order=order,
+ boundary_tag_to_face=boundary_tag_to_face,
+ group_factory=group_factory)
# }}}
diff --git a/meshmode/mesh/processing.py b/meshmode/mesh/processing.py
index 36e2505e8617887e95d37cfaf79b0eb07eb413ee..70c17461bc5050e74b2e640bd22b2f034eb1c60a 100644
--- a/meshmode/mesh/processing.py
+++ b/meshmode/mesh/processing.py
@@ -37,6 +37,7 @@ __doc__ = """
.. autofunction:: perform_flips
.. autofunction:: find_bounding_box
.. autofunction:: merge_disjoint_meshes
+.. autofunction:: split_mesh_groups
.. autofunction:: map_mesh
.. autofunction:: affine_map
"""
@@ -588,6 +589,64 @@ def merge_disjoint_meshes(meshes, skip_tests=False, single_group=False):
# }}}
+# {{{ split meshes
+
+def split_mesh_groups(mesh, element_flags, return_subgroup_mapping=False):
+ """Split all the groups in *mesh* according to the values of
+ *element_flags*. The element flags are expected to be integers
+ defining, for each group, how the elements are to be split into
+ subgroups. For example, a single-group mesh with flags::
+
+ element_flags = [0, 0, 0, 42, 42, 42, 0, 0, 0, 41, 41, 41]
+
+ will create three subgroups. The integer flags need not be increasing
+ or contiguous and can repeat across different groups (i.e. they are
+ group-local).
+
+ :arg element_flags: a :class:`numpy.ndarray` with
+ :attr:`~meshmode.mesh.Mesh.nelements` entries
+ indicating how the elements in a group are to be split.
+
+ :returns: a :class:`~meshmode.mesh.Mesh` where each group has been split
+ according to flags in *element_flags*. If *return_subgroup_mapping*
+ is *True*, it also returns a mapping of
+ ``(group_index, subgroup) -> new_group_index``.
+
+ """
+ assert element_flags.shape == (mesh.nelements,)
+
+ new_groups = []
+ subgroup_to_group_map = {}
+
+ for igrp, grp in enumerate(mesh.groups):
+ grp_flags = element_flags[
+ grp.element_nr_base:grp.element_nr_base + grp.nelements]
+ unique_grp_flags = np.unique(grp_flags)
+
+ for flag in unique_grp_flags:
+ subgroup_to_group_map[igrp, flag] = len(new_groups)
+
+ # NOTE: making copies to maintain contiguity of the arrays
+ mask = grp_flags == flag
+ new_groups.append(grp.copy(
+ vertex_indices=grp.vertex_indices[mask, :].copy(),
+ nodes=grp.nodes[:, mask, :].copy()
+ ))
+
+ from meshmode.mesh import Mesh
+ mesh = Mesh(
+ vertices=mesh.vertices,
+ groups=new_groups,
+ is_conforming=mesh.is_conforming)
+
+ if return_subgroup_mapping:
+ return mesh, subgroup_to_group_map
+ else:
+ return mesh
+
+# }}}
+
+
# {{{ map
def map_mesh(mesh, f): # noqa
diff --git a/meshmode/mesh/refinement/no_adjacency.py b/meshmode/mesh/refinement/no_adjacency.py
index f37f66215ded63f530db8d5a2cb465d38771b5db..87f843d3ac052ebbb39b20d85deeb76445b2c5db 100644
--- a/meshmode/mesh/refinement/no_adjacency.py
+++ b/meshmode/mesh/refinement/no_adjacency.py
@@ -238,6 +238,8 @@ class RefinerWithoutAdjacency(object):
global_midpoint = inew_vertex
additional_vertices.append(
midpoints[old_iel][:, imidpoint])
+ self.global_vertex_pair_to_midpoint[
+ global_v1, global_v2] = global_midpoint
inew_vertex += 1
refining_vertices[iref_midpoint] = global_midpoint
diff --git a/meshmode/version.py b/meshmode/version.py
index d26fbc2f9341a880b1119e7e6079bd51e59e11b9..5cea8857d9550183c2caabd0948032faac48e0ae 100644
--- a/meshmode/version.py
+++ b/meshmode/version.py
@@ -1,2 +1,2 @@
-VERSION = (2016, 1)
+VERSION = (2020, 1)
VERSION_TEXT = ".".join(str(i) for i in VERSION)
diff --git a/setup.py b/setup.py
index 688e9f92c8dbf573b4b5702da2cd50b420ca08d1..55a523ca5cf961abd0c8bfb9b70c35f0898b0b5e 100644
--- a/setup.py
+++ b/setup.py
@@ -27,10 +27,6 @@ def main():
'Natural Language :: English',
'Programming Language :: Python',
'Programming Language :: Python :: 3',
- 'Programming Language :: Python :: 2.6',
- 'Programming Language :: Python :: 2.7',
- 'Programming Language :: Python :: 3.2',
- 'Programming Language :: Python :: 3.3',
'Topic :: Scientific/Engineering',
'Topic :: Scientific/Engineering :: Information Analysis',
'Topic :: Scientific/Engineering :: Mathematics',
@@ -40,6 +36,7 @@ def main():
],
packages=find_packages(),
+ python_requires="~=3.6",
install_requires=[
"numpy",
"modepy",
diff --git a/test/test_meshmode.py b/test/test_meshmode.py
index 77a25cd5dd4de07e1207ca137ec5b1510d268a84..047362d20c1b297e22eed758fbe633bee01a82cc 100644
--- a/test/test_meshmode.py
+++ b/test/test_meshmode.py
@@ -33,12 +33,13 @@ from pyopencl.tools import ( # noqa
pytest_generate_tests_for_pyopencl
as pytest_generate_tests)
+from meshmode.mesh import SimplexElementGroup, TensorProductElementGroup
from meshmode.discretization.poly_element import (
InterpolatoryQuadratureSimplexGroupFactory,
PolynomialWarpAndBlendGroupFactory,
PolynomialEquidistantSimplexGroupFactory,
)
-from meshmode.mesh import BTAG_ALL
+from meshmode.mesh import Mesh, BTAG_ALL
from meshmode.discretization.connection import \
FACE_RESTR_ALL, FACE_RESTR_INTERIOR
import meshmode.mesh.generation as mgen
@@ -123,6 +124,84 @@ def test_boundary_tags():
# }}}
+# {{{ test custom boundary tags on box mesh
+
+@pytest.mark.parametrize(("dim", "nelem"), [
+ (1, 20),
+ (2, 20),
+ (3, 10),
+ ])
+@pytest.mark.parametrize("group_factory", [
+ SimplexElementGroup,
+
+ # FIXME: Not implemented: TPE.face_vertex_indices
+ # TensorProductElementGroup
+ ])
+def test_box_boundary_tags(dim, nelem, group_factory):
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ from meshmode.mesh import is_boundary_tag_empty
+ from meshmode.mesh import check_bc_coverage
+ if dim == 1:
+ a = (0,)
+ b = (1,)
+ n = (nelem,)
+ btag_to_face = {"btag_test_1": ["+x"],
+ "btag_test_2": ["-x"]}
+ elif dim == 2:
+ a = (0, -1)
+ b = (1, 1)
+ n = (nelem, nelem)
+ btag_to_face = {"btag_test_1": ["+x", "-y"],
+ "btag_test_2": ["+y", "-x"]}
+ elif dim == 3:
+ a = (0, -1, -1)
+ b = (1, 1, 1)
+ n = (nelem, nelem, nelem)
+ btag_to_face = {"btag_test_1": ["+x", "-y", "-z"],
+ "btag_test_2": ["+y", "-x", "+z"]}
+ mesh = generate_regular_rect_mesh(a=a, b=b,
+ n=n, order=3,
+ boundary_tag_to_face=btag_to_face,
+ group_factory=group_factory)
+ # correct answer
+ if dim == 1:
+ num_on_bdy = 1
+ else:
+ num_on_bdy = dim*(dim-1)*(nelem-1)**(dim-1)
+
+ assert not is_boundary_tag_empty(mesh, "btag_test_1")
+ assert not is_boundary_tag_empty(mesh, "btag_test_2")
+ check_bc_coverage(mesh, ['btag_test_1', 'btag_test_2'])
+
+ # check how many elements are marked on each boundary
+ num_marked_bdy_1 = 0
+ num_marked_bdy_2 = 0
+ btag_1_bit = mesh.boundary_tag_bit("btag_test_1")
+ btag_2_bit = mesh.boundary_tag_bit("btag_test_2")
+ for igrp in range(len(mesh.groups)):
+ bdry_fagrp = mesh.facial_adjacency_groups[igrp].get(None, None)
+
+ if bdry_fagrp is None:
+ continue
+
+ for i, nbrs in enumerate(bdry_fagrp.neighbors):
+ if (-nbrs) & btag_1_bit:
+ num_marked_bdy_1 += 1
+ if (-nbrs) & btag_2_bit:
+ num_marked_bdy_2 += 1
+
+ # raise errors if wrong number of elements marked
+ if num_marked_bdy_1 != num_on_bdy:
+ raise ValueError("%i marked on custom boundary 1, should be %i" %
+ (num_marked_bdy_1, num_on_bdy))
+ if num_marked_bdy_2 != num_on_bdy:
+ raise ValueError("%i marked on custom boundary 2, should be %i" %
+ (num_marked_bdy_2, num_on_bdy))
+
+
+# }}}
+
+
# {{{ convergence of boundary interpolation
@pytest.mark.parametrize("group_factory", [
@@ -527,7 +606,6 @@ def test_3d_orientation(ctx_factory, what, mesh_gen_func, visualize=False):
def test_merge_and_map(ctx_factory, visualize=False):
from meshmode.mesh.io import generate_gmsh, FileSource
from meshmode.mesh.generation import generate_box_mesh
- from meshmode.mesh import TensorProductElementGroup
from meshmode.discretization.poly_element import (
PolynomialWarpAndBlendGroupFactory,
LegendreGaussLobattoTensorProductGroupFactory)
@@ -1009,7 +1087,6 @@ def no_test_quad_mesh_3d():
def test_quad_single_element():
from meshmode.mesh.generation import make_group_from_vertices
- from meshmode.mesh import Mesh, TensorProductElementGroup
vertices = np.array([
[0.91, 1.10],
@@ -1037,7 +1114,6 @@ def test_quad_single_element():
def test_quad_multi_element():
from meshmode.mesh.generation import generate_box_mesh
- from meshmode.mesh import TensorProductElementGroup
mesh = generate_box_mesh(
(
np.linspace(3, 8, 4),
@@ -1272,6 +1348,89 @@ def test_is_affine_group_check(mesh_name):
assert all(grp.is_affine for grp in mesh.groups) == is_affine
+@pytest.mark.parametrize("ambient_dim", [1, 2, 3])
+def test_mesh_multiple_groups(ctx_factory, ambient_dim, visualize=False):
+ ctx = ctx_factory()
+ queue = cl.CommandQueue(ctx)
+
+ order = 4
+
+ from meshmode.mesh.generation import generate_regular_rect_mesh
+ mesh = generate_regular_rect_mesh(
+ a=(-0.5,)*ambient_dim, b=(0.5,)*ambient_dim,
+ n=(8,)*ambient_dim, order=order)
+ assert len(mesh.groups) == 1
+
+ from meshmode.mesh.processing import split_mesh_groups
+ element_flags = np.any(
+ mesh.vertices[0, mesh.groups[0].vertex_indices] < 0.0,
+ axis=1).astype(np.int)
+ mesh = split_mesh_groups(mesh, element_flags)
+
+ assert len(mesh.groups) == 2
+ assert mesh.facial_adjacency_groups
+ assert mesh.nodal_adjacency
+
+ if visualize and ambient_dim == 2:
+ from meshmode.mesh.visualization import draw_2d_mesh
+ draw_2d_mesh(mesh,
+ draw_vertex_numbers=False,
+ draw_element_numbers=True,
+ draw_face_numbers=False,
+ set_bounding_box=True)
+
+ import matplotlib.pyplot as plt
+ plt.savefig("test_mesh_multiple_groups_2d_elements.png", dpi=300)
+
+ from meshmode.discretization import Discretization
+ from meshmode.discretization.poly_element import \
+ PolynomialWarpAndBlendGroupFactory as GroupFactory
+ discr = Discretization(ctx, mesh, GroupFactory(order))
+
+ if visualize:
+ group_id = discr.empty(queue, dtype=np.int)
+ for igrp, grp in enumerate(discr.groups):
+ group_id_view = grp.view(group_id)
+ group_id_view.fill(igrp)
+
+ from meshmode.discretization.visualization import make_visualizer
+ vis = make_visualizer(queue, discr, vis_order=order)
+ vis.write_vtk_file("test_mesh_multiple_groups.vtu", [
+ ("group_id", group_id)
+ ], overwrite=True)
+
+ # check face restrictions
+ from meshmode.discretization.connection import (
+ make_face_restriction,
+ make_face_to_all_faces_embedding,
+ make_opposite_face_connection,
+ check_connection)
+ for boundary_tag in [BTAG_ALL, FACE_RESTR_INTERIOR, FACE_RESTR_ALL]:
+ conn = make_face_restriction(discr, GroupFactory(order),
+ boundary_tag=boundary_tag,
+ per_face_groups=False)
+ check_connection(conn)
+
+ bdry_f = conn.to_discr.empty(queue)
+ bdry_f.fill(1.0)
+
+ if boundary_tag == FACE_RESTR_INTERIOR:
+ opposite = make_opposite_face_connection(conn)
+ check_connection(opposite)
+
+ op_bdry_f = opposite(queue, bdry_f)
+ error = abs(cl.array.sum(bdry_f - op_bdry_f).get(queue))
+ assert error < 1.0e-11, error
+
+ if boundary_tag == FACE_RESTR_ALL:
+ embedding = make_face_to_all_faces_embedding(conn, conn.to_discr)
+ check_connection(embedding)
+
+ em_bdry_f = embedding(queue, bdry_f)
+ error = abs(cl.array.sum(bdry_f - em_bdry_f).get(queue))
+ assert error < 1.0e-11, error
+
+
if __name__ == "__main__":
import sys
if len(sys.argv) > 1:
diff --git a/test/test_refinement.py b/test/test_refinement.py
index 5e10a9462ae7b1aa195227a58950ec69482775ff..2bc29e89266f4e7a5ae824335d719cb259412ea6 100644
--- a/test/test_refinement.py
+++ b/test/test_refinement.py
@@ -294,7 +294,7 @@ def test_refinement_connection(
@pytest.mark.parametrize("with_adjacency", [True, False])
-def test_uniform_refinement(ctx_factory, with_adjacency):
+def test_uniform_refinement(with_adjacency):
make_mesh = partial(generate_box_mesh, (
np.linspace(0.0, 1.0, 2),
np.linspace(0.0, 1.0, 3),
@@ -305,6 +305,18 @@ def test_uniform_refinement(ctx_factory, with_adjacency):
mesh = refine_uniformly(mesh, 1, with_adjacency=with_adjacency)
+@pytest.mark.parametrize("refinement_rounds", [0, 1, 2])
+def test_conformity_of_uniform_mesh(refinement_rounds):
+ from meshmode.mesh.generation import generate_icosphere
+ mesh = generate_icosphere(r=1.0, order=4,
+ uniform_refinement_rounds=refinement_rounds)
+
+ assert mesh.is_conforming
+
+ from meshmode.mesh import is_boundary_tag_empty, BTAG_ALL
+ assert is_boundary_tag_empty(mesh, BTAG_ALL)
+
+
if __name__ == "__main__":
import sys
if len(sys.argv) > 1: