From 8406098dfbbc58077276ae09d80436b101633ac7 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 5 Jun 2016 13:03:37 -0500 Subject: [PATCH 01/23] Start working on quad meshes --- meshmode/mesh/io.py | 14 ++++++++++---- test/test_meshmode.py | 23 +++++++++++++++++++++++ 2 files changed, 33 insertions(+), 4 deletions(-) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index e27f8de3..4a97c72e 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -27,13 +27,15 @@ from six.moves import range, zip import numpy as np from meshpy.gmsh_reader import ( # noqa - GmshMeshReceiverBase, FileSource, LiteralSource) + GmshMeshReceiverBase, ScriptSource, FileSource, LiteralSource, + ScriptWithFilesSource) __doc__ = """ +.. autoclass:: ScriptSource .. autoclass:: FileSource -.. autoclass:: LiteralSource +.. autoclass:: ScriptWithFilesSource .. autofunction:: read_gmsh .. autofunction:: generate_gmsh @@ -96,6 +98,9 @@ class GmshMeshReceiver(GmshMeshReceiverBase): for el_type in self.element_types: el_type_hist[el_type] = el_type_hist.get(el_type, 0) + 1 + if not el_type_hist: + raise RuntimeError("empty mesh in gmsh input") + groups = self.groups = [] ambient_dim = self.points.shape[-1] @@ -195,8 +200,9 @@ def read_gmsh(filename, force_ambient_dim=None): return recv.get_mesh() -def generate_gmsh(source, dimensions, order=None, other_options=[], - extension="geo", gmsh_executable="gmsh", force_ambient_dim=None): +def generate_gmsh(source, dimensions=None, order=None, other_options=[], + extension="geo", gmsh_executable="gmsh", force_ambient_dim=None, + output_file_name="output.msh"): """Run :command:`gmsh` on the input given by *source*, and return a :class:`meshmode.mesh.Mesh` based on the result. diff --git a/test/test_meshmode.py b/test/test_meshmode.py index 5eba7feb..b3ca2022 100644 --- a/test/test_meshmode.py +++ b/test/test_meshmode.py @@ -802,6 +802,29 @@ def test_nd_quad_submesh(dims): # }}} +# {{{ test_quad_mesh + +def test_quad_mesh(): + from meshmode.mesh.io import generate_gmsh, ScriptWithFilesSource + print("BEGIN GEN") + mesh = generate_gmsh( + ScriptWithFilesSource( + """ + Merge "blob-2d.step"; + Mesh.CharacteristicLengthMax = 0.05; + Recombine Surface "*" = 0.0001; + Mesh 2; + Save "output.msh"; + """, + ["blob-2d.step"]), + force_ambient_dim=2, + ) + print("END GEN") + print(mesh.nelements) + +# }}} + + if __name__ == "__main__": import sys if len(sys.argv) > 1: -- GitLab From 44a6e162e7bf287763b71c4e4bf974dacfb22615 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Sun, 12 Jun 2016 00:01:37 -0500 Subject: [PATCH 02/23] Bump Py3 CI tests to Py 3.5 --- .gitlab-ci.yml | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 609ef142..bc7599d4 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -22,15 +22,15 @@ Python 2.7 POCL: - pocl except: - tags -Python 3.4 POCL: +Python 3.5 POCL: script: - - export PY_EXE=python3.4 + - export PY_EXE=python3.5 - export PYOPENCL_TEST=portable - export EXTRA_INSTALL="numpy mako" - 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: - - python3.4 + - python3.5 - pocl except: - tags -- GitLab From 0822ba0c2a8954859040a3719c0021eac0381b03 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Tue, 28 Jun 2016 20:13:52 -0500 Subject: [PATCH 03/23] Initial support for quad meshes, also in Gmsh I/O --- meshmode/mesh/__init__.py | 36 ++++++++++++++++++++++++++++++++++ meshmode/mesh/io.py | 41 ++++++++++++++++++++++++++------------- 2 files changed, 64 insertions(+), 13 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 5c2b5c1c..4e2f1f0a 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -292,6 +292,42 @@ class SimplexElementGroup(MeshElementGroup): from modepy.tools import unit_vertices return unit_vertices(self.dim) + +class TensorProductElementGroup(MeshElementGroup): + def __init__(self, order, vertex_indices, nodes, + element_nr_base=None, node_nr_base=None, + unit_nodes=None): + """ + :arg order: the mamximum total degree used for interpolation. + :arg nodes: ``[ambient_dim, nelements, nunit_nodes]`` + The nodes are assumed to be mapped versions of *unit_nodes*. + :arg unit_nodes: ``[dim, nunit_nodes]`` + The unit nodes of which *nodes* is a mapped + version. + + Do not supply *element_nr_base* and *node_nr_base*, they will be + automatically assigned. + """ + + if not issubclass(vertex_indices.dtype.type, np.integer): + raise TypeError("vertex_indices must be integral") + + dims = unit_nodes.shape[0] + + if vertex_indices.shape[-1] != 2**dims: + raise ValueError("vertex_indices has wrong number of vertices per " + "element. expected: %d, got: %d" % (dims+1, + vertex_indices.shape[-1])) + + super(TensorProductElementGroup, self).__init__(order, vertex_indices, nodes, + element_nr_base, node_nr_base, unit_nodes) + + def face_vertex_indices(self): + raise NotImplementedError() + + def vertex_unit_coordinates(self): + raise NotImplementedError() + # }}} diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index 4a97c72e..a5bbe204 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -28,7 +28,9 @@ import numpy as np from meshpy.gmsh_reader import ( # noqa GmshMeshReceiverBase, ScriptSource, FileSource, LiteralSource, - ScriptWithFilesSource) + ScriptWithFilesSource, + GmshSimplexElementBase, + GmshTensorProductElementBase) __doc__ = """ @@ -131,7 +133,8 @@ class GmshMeshReceiver(GmshMeshReceiverBase): # }}} - from meshmode.mesh import SimplexElementGroup, Mesh + from meshmode.mesh import (Mesh, + SimplexElementGroup, TensorProductElementGroup) for group_el_type, ngroup_elements in six.iteritems(el_type_hist): if group_el_type.dimensions != mesh_bulk_dim: @@ -159,21 +162,33 @@ class GmshMeshReceiver(GmshMeshReceiverBase): unit_nodes = (np.array(group_el_type.lexicographic_node_tuples(), dtype=np.float64).T/group_el_type.order)*2 - 1 - group = SimplexElementGroup( - group_el_type.order, - vertex_indices, - nodes, - unit_nodes=unit_nodes - ) + if isinstance(group_el_type, GmshSimplexElementBase): + group = SimplexElementGroup( + group_el_type.order, + vertex_indices, + nodes, + unit_nodes=unit_nodes + ) + + if group.dim == 2: + from meshmode.mesh.processing import flip_simplex_element_group + group = flip_simplex_element_group(vertices, group, + np.ones(ngroup_elements, np.bool)) + + elif isinstance(group_el_type, GmshTensorProductElementBase): + group = TensorProductElementGroup( + group_el_type.order, + vertex_indices, + nodes, + unit_nodes=unit_nodes + ) + else: + raise NotImplementedError("gmsh element type: %s" + % type(group_el_type).__name__) # Gmsh seems to produce elements in the opposite orientation # of what we like. Flip them all. - if group.dim == 2: - from meshmode.mesh.processing import flip_simplex_element_group - group = flip_simplex_element_group(vertices, group, - np.ones(ngroup_elements, np.bool)) - groups.append(group) return Mesh( -- GitLab From e434f44fca1aed74560f3b16112b0276b373470e Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 29 Jun 2016 12:28:00 -0500 Subject: [PATCH 04/23] Requirements: switch to https, use meshpy from git --- requirements.txt | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/requirements.txt b/requirements.txt index 8e99cd7f..0a0cff88 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,10 +1,11 @@ numpy -git+git://github.com/inducer/modepy -git+git://github.com/pyopencl/pyopencl -git+git://github.com/inducer/islpy -git+git://github.com/inducer/loopy +https+git://github.com/inducer/meshpy.git +https+git://github.com/inducer/modepy.git +https+git://github.com/pyopencl/pyopencl.git +https+git://github.com/inducer/islpy.git +https+git://github.com/inducer/loopy.git # required by pytential, which is in turn needed for some tests -git+git://github.com/inducer/boxtree -git+git://github.com/inducer/sumpy -git+git://github.com/inducer/pytential +https+git://github.com/inducer/boxtree.git +https+git://github.com/inducer/sumpy.git +https+git://github.com/inducer/pytential.git -- GitLab From d4513a0c5005a756ee5e88d90bb806b91c6a6dcb Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 29 Jun 2016 12:29:46 -0500 Subject: [PATCH 05/23] Fix requirements git URLs --- requirements.txt | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/requirements.txt b/requirements.txt index 0a0cff88..6ed102e6 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,11 +1,11 @@ numpy -https+git://github.com/inducer/meshpy.git -https+git://github.com/inducer/modepy.git -https+git://github.com/pyopencl/pyopencl.git -https+git://github.com/inducer/islpy.git -https+git://github.com/inducer/loopy.git +git+https://github.com/inducer/meshpy.git +git+https://github.com/inducer/modepy.git +git+https://github.com/pyopencl/pyopencl.git +git+https://github.com/inducer/islpy.git +git+https://github.com/inducer/loopy.git # required by pytential, which is in turn needed for some tests -https+git://github.com/inducer/boxtree.git -https+git://github.com/inducer/sumpy.git -https+git://github.com/inducer/pytential.git +git+https://github.com/inducer/boxtree.git +git+https://github.com/inducer/sumpy.git +git+https://github.com/inducer/pytential.git -- GitLab From 2563ab8092b6f788fa3384a5728b07b5c6d7d424 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Thu, 30 Jun 2016 10:39:53 -0500 Subject: [PATCH 06/23] Mini loopy modernization --- meshmode/discretization/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/discretization/__init__.py b/meshmode/discretization/__init__.py index e1ed150f..45fdc22e 100644 --- a/meshmode/discretization/__init__.py +++ b/meshmode/discretization/__init__.py @@ -311,7 +311,7 @@ class Discretization(object): knl = lp.split_iname(knl, "i", 16, inner_tag="l.0") knl = lp.tag_inames(knl, dict(k="g.0")) - knl = lp.tag_data_axes(knl, "result", + knl = lp.tag_array_axes(knl, "result", "stride:auto,stride:auto,stride:auto") return knl -- GitLab From 09927b549a156a70bf9a7462bb9aacf331e59901 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 8 Aug 2016 11:47:06 -0500 Subject: [PATCH 07/23] Make refiner-generated data structure match docs (fix attr typo) --- meshmode/mesh/refinement/__init__.py | 2 +- meshmode/mesh/refinement/utils.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index 2154d451..ac2f9833 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -845,7 +845,7 @@ class Refiner(object): assert neighbors_starts[-1] == len(neighbors) from meshmode.mesh import NodalAdjacency - return NodalAdjacency(neighbor_starts=neighbors_starts, neighbors=neighbors) + return NodalAdjacency(neighbors_starts=neighbors_starts, neighbors=neighbors) # }}} diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index 9a78c31d..0d0476f3 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -64,7 +64,7 @@ def check_nodal_adj_against_geometry(mesh, tol=1e-12): for igrp, grp in enumerate(mesh.groups): for iel_grp in range(grp.nelements): iel_g = group_and_iel_to_global_iel(igrp, iel_grp) - nb_starts = nadj.neighbor_starts + nb_starts = nadj.neighbors_starts for nb_iel_g in nadj.neighbors[nb_starts[iel_g]:nb_starts[iel_g+1]]: connected_to_element_connectivity[iel_g].add(nb_iel_g) -- GitLab From 1987619bf10c14b53664af8640d234721176e07c Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Fri, 29 Jul 2016 16:23:45 -0500 Subject: [PATCH 08/23] make_refinement_connection(): Improve error behavior. --- meshmode/discretization/connection/refinement.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 75604ec0..4477a9bf 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -161,7 +161,9 @@ def make_refinement_connection(refiner, coarse_discr, group_factory): coarse_mesh = refiner.get_previous_mesh() fine_mesh = refiner.last_mesh - assert coarse_discr.mesh is coarse_mesh + if coarse_discr.mesh != coarse_mesh: + raise ValueError( + "course_discr must live on the same mesh given to the refiner!") from meshmode.discretization import Discretization fine_discr = Discretization( -- GitLab From 7d76f89337a516713eb039ddd454b825a4578e15 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 8 Aug 2016 11:47:49 -0500 Subject: [PATCH 09/23] Bump version --- meshmode/version.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/meshmode/version.py b/meshmode/version.py index 1da70a1e..d26fbc2f 100644 --- a/meshmode/version.py +++ b/meshmode/version.py @@ -1,2 +1,2 @@ -VERSION = (2015, 1) +VERSION = (2016, 1) VERSION_TEXT = ".".join(str(i) for i in VERSION) -- GitLab From 32fd5f8689d7bed5e76871eb4fc071a1ce8f2e0a Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 8 Aug 2016 11:48:40 -0500 Subject: [PATCH 10/23] Add to_json mesh output --- meshmode/mesh/io.py | 54 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 54 insertions(+) diff --git a/meshmode/mesh/io.py b/meshmode/mesh/io.py index a5bbe204..f3cdaea8 100644 --- a/meshmode/mesh/io.py +++ b/meshmode/mesh/io.py @@ -43,6 +43,7 @@ __doc__ = """ .. autofunction:: generate_gmsh .. autofunction:: from_meshpy .. autofunction:: from_vertices_and_simplices +.. autofunction:: to_json """ @@ -314,4 +315,57 @@ def from_vertices_and_simplices(vertices, simplices, order=1, fix_orientation=Fa # }}} +# {{{ to_json + +def to_json(mesh): + """Return a JSON-able Python data structure for *mesh*. The structure directly + reflects the :class:`Mesh` data structure.""" + + def btag_to_json(btag): + if isinstance(btag, str): + return btag + else: + return btag.__name__ + + def group_to_json(group): + return { + "type": type(group).__name__, + "order": group.order, + "vertex_indices": group.vertex_indices.tolist(), + "nodes": group.nodes.tolist(), + "unit_nodes": group.unit_nodes.tolist(), + "element_nr_base": group.element_nr_base, + "node_nr_base": group.node_nr_base, + "dim": group.dim, + } + + from meshmode import DataUnavailable + + def nodal_adjacency_to_json(mesh): + try: + na = mesh.nodal_adjacency + except DataUnavailable: + return None + + return { + "neighbors_starts": na.neighbors_starts.tolist(), + "neighbors": na.neighbors.tolist(), + } + + return { + "version": 0, + "vertices": mesh.vertices.tolist(), + "groups": [group_to_json(group) for group in mesh.groups], + "nodal_adjacency": nodal_adjacency_to_json(mesh), + # not yet implemented + "facial_adjacency_groups": None, + "boundary_tags": [btag_to_json(btag) for btag in mesh.boundary_tags], + "btag_to_index": dict( + (btag_to_json(btag), value) + for btag, value in six.iteritems(mesh.btag_to_index)), + } + +# }}} + + # vim: foldmethod=marker -- GitLab From e0aab6ff8c67477a6aa73810a0d58384470cee03 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 8 Aug 2016 11:49:01 -0500 Subject: [PATCH 11/23] Add misc doc section --- doc/index.rst | 1 + doc/misc.rst | 79 +++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 80 insertions(+) create mode 100644 doc/misc.rst diff --git a/doc/index.rst b/doc/index.rst index f652d1ba..4bde9b27 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -8,6 +8,7 @@ Contents: mesh discretization + misc Indices and tables ================== diff --git a/doc/misc.rst b/doc/misc.rst new file mode 100644 index 00000000..488262ac --- /dev/null +++ b/doc/misc.rst @@ -0,0 +1,79 @@ +.. _installation: + +Installation +============ + +This command should install :mod:`meshmode`:: + + pip install meshmode + +(Note the extra "."!) + +You may need to run this with :command:`sudo`. +If you don't already have `pip `_, +run this beforehand:: + + curl -O https://raw.github.com/pypa/pip/master/contrib/get-pip.py + python get-pip.py + +For a more manual installation, `download the source +`_, unpack it, and say:: + + python setup.py install + +You may also clone its git repository:: + + git clone --recursive git://github.com/inducer/meshmode + git clone --recursive http://git.tiker.net/trees/meshmode.git + +User-visible Changes +==================== + +Version 2016.1 +-------------- +.. note:: + + This version is currently under development. You can get snapshots from + meshmode's `git repository `_ + +.. _license: + +Licensing +========= + +:mod:`meshmode` is licensed to you under the MIT/X Consortium license: + +Copyright (c) 2014-16 Andreas Klöckner and Contributors. + +Permission is hereby granted, free of charge, to any person +obtaining a copy of this software and associated documentation +files (the "Software"), to deal in the Software without +restriction, including without limitation the rights to use, +copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the +Software is furnished to do so, subject to the following +conditions: + +The above copyright notice and this permission notice shall be +included in all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES +OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND +NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT +HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, +WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING +FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR +OTHER DEALINGS IN THE SOFTWARE. + +Acknowledgments +=============== + +Andreas Klöckner's work on :mod:`meshmode` was supported in part by + +* US Navy ONR grant number N00014-14-1-0117 +* the US National Science Foundation under grant numbers DMS-1418961 and CCF-1524433. + +AK also gratefully acknowledges a hardware gift from Nvidia Corporation. The +views and opinions expressed herein do not necessarily reflect those of the +funding agencies. -- GitLab From 7d2e14422121e65ec9028dd9f81a7429c8526f33 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 8 Aug 2016 11:49:15 -0500 Subject: [PATCH 12/23] Add doc autobuild --- .gitlab-ci.yml | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index bc7599d4..7c2f2c97 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -34,3 +34,12 @@ Python 3.5 POCL: - pocl except: - tags +Documentation: + script: + - EXTRA_INSTALL="numpy" + - curl -L -O -k https://gitlab.tiker.net/inducer/ci-support/raw/master/build-docs.sh + - ". ./build-docs.sh" + tags: + - python3.5 + only: + - master -- GitLab From dde79a7bada6156c8018a3c80963d9de4b623f48 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Mon, 8 Aug 2016 12:20:36 -0500 Subject: [PATCH 13/23] Fix more neighbor_starts -> neighbors_starts confusion --- meshmode/discretization/visualization.py | 2 +- meshmode/mesh/__init__.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/meshmode/discretization/visualization.py b/meshmode/discretization/visualization.py index 956e3015..d68ec0c5 100644 --- a/meshmode/discretization/visualization.py +++ b/meshmode/discretization/visualization.py @@ -291,7 +291,7 @@ def write_nodal_adjacency_vtk_file(file_name, mesh, compressor=None,): nconnections = len(adj.neighbors) connections = np.empty((nconnections, 2), dtype=np.int32) - nb_starts = adj.neighbor_starts + nb_starts = adj.neighbors_starts for iel in range(mesh.nelements): connections[nb_starts[iel]:nb_starts[iel+1], 0] = iel diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 4e2f1f0a..38286064 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -581,10 +581,10 @@ class Mesh(Record): assert g.vertex_indices.dtype == self.vertex_id_dtype if nodal_adjacency: - assert nodal_adjacency.neighbor_starts.shape == (self.nelements+1,) + assert nodal_adjacency.neighbors_starts.shape == (self.nelements+1,) assert len(nodal_adjacency.neighbors.shape) == 1 - assert nodal_adjacency.neighbor_starts.dtype == self.element_id_dtype + assert nodal_adjacency.neighbors_starts.dtype == self.element_id_dtype assert nodal_adjacency.neighbors.dtype == self.element_id_dtype if facial_adjacency_groups: -- GitLab From e1c1ee74bb44a55e9d463567f430ad890d4e6fff Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 10 Aug 2016 18:16:56 -0500 Subject: [PATCH 14/23] Ignore pytest's cache directory --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index feac29bf..e090b1d8 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,5 @@ distribute*egg distribute*tar.gz a.out *.vtu + +.cache -- GitLab From 8266509059911ba7dacd87ea484110fced184107 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 10 Aug 2016 18:16:59 -0500 Subject: [PATCH 15/23] Don't show progress during doc upload --- doc/upload-docs.sh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/upload-docs.sh b/doc/upload-docs.sh index 25c97f52..2c219a8c 100755 --- a/doc/upload-docs.sh +++ b/doc/upload-docs.sh @@ -1,3 +1,3 @@ #! /bin/sh -rsync --progress --verbose --archive --delete _build/html/* doc-upload:doc/meshmode +rsync --verbose --archive --delete _build/html/* doc-upload:doc/meshmode -- GitLab From d5cd39d6fe990585f09e6b8a5e42da7d250b78c0 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Wed, 10 Aug 2016 18:17:09 -0500 Subject: [PATCH 16/23] PEP8 cleanups --- meshmode/mesh/visualization.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/meshmode/mesh/visualization.py b/meshmode/mesh/visualization.py index caa9c1be..601086c1 100644 --- a/meshmode/mesh/visualization.py +++ b/meshmode/mesh/visualization.py @@ -1,6 +1,4 @@ -from __future__ import division -from __future__ import absolute_import -from six.moves import range +from __future__ import division, absolute_import __copyright__ = "Copyright (C) 2014 Andreas Kloeckner" @@ -25,6 +23,7 @@ THE SOFTWARE. """ import numpy as np +from six.moves import range # {{{ draw_2d_mesh -- GitLab From 4623caeee603eee0cd0582da640f5f6c3eb68ef8 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 12 Aug 2016 14:25:00 -0500 Subject: [PATCH 17/23] Mesh: Add node_vertex_consistency_tolerance --- meshmode/mesh/__init__.py | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/meshmode/mesh/__init__.py b/meshmode/mesh/__init__.py index 38286064..da0fdfe2 100644 --- a/meshmode/mesh/__init__.py +++ b/meshmode/mesh/__init__.py @@ -490,6 +490,7 @@ class Mesh(Record): face_id_dtype = np.int8 def __init__(self, vertices, groups, skip_tests=False, + node_vertex_consistency_tolerance=None, nodal_adjacency=False, facial_adjacency_groups=False, boundary_tags=None, @@ -500,6 +501,9 @@ class Mesh(Record): :arg skip_tests: Skip mesh tests, in case you want to load a broken mesh anyhow and then fix it inside of this data structure. + :arg node_vertex_consistency_tolerance: If *False*, do not check + for consistency between vertex and nodal data. If *None*, use + the (small, near FP-epsilon) tolerance. :arg nodal_adjacency: One of three options: *None*, in which case this information will be deduced from vertex adjacency. *False*, in which case @@ -576,7 +580,8 @@ class Mesh(Record): ) if not skip_tests: - assert _test_node_vertex_consistency(self) + assert _test_node_vertex_consistency( + self, node_vertex_consistency_tolerance) for g in self.groups: assert g.vertex_indices.dtype == self.vertex_id_dtype @@ -691,7 +696,7 @@ class Mesh(Record): # {{{ node-vertex consistency test -def _test_node_vertex_consistency_simplex(mesh, mgrp): +def _test_node_vertex_consistency_simplex(mesh, mgrp, tol): if mgrp.nelements == 0: return True @@ -710,7 +715,8 @@ def _test_node_vertex_consistency_simplex(mesh, mgrp): np.sum((map_vertices - grp_vertices)**2, axis=0), axis=-1)) - tol = 1e3 * np.finfo(per_element_vertex_errors.dtype).eps + if tol is None: + tol = 1e3 * np.finfo(per_element_vertex_errors.dtype).eps from meshmode.mesh.processing import find_bounding_box @@ -723,14 +729,17 @@ def _test_node_vertex_consistency_simplex(mesh, mgrp): return True -def _test_node_vertex_consistency(mesh): +def _test_node_vertex_consistency(mesh, tol): """Ensure that order of by-index vertices matches that of mapped unit vertices. """ + if tol is False: + return + for mgrp in mesh.groups: if isinstance(mgrp, SimplexElementGroup): - assert _test_node_vertex_consistency_simplex(mesh, mgrp) + assert _test_node_vertex_consistency_simplex(mesh, mgrp, tol) else: from warnings import warn warn("not implemented: node-vertex consistency check for '%s'" -- GitLab From d9f38e45903927c5f90ba566067a795a818404c2 Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 12 Aug 2016 14:25:27 -0500 Subject: [PATCH 18/23] make_curve_mesh: Allow specifying unit nodes and node_vertex_consistency_tolerance --- meshmode/mesh/generation.py | 31 +++++++++++++++++++++++-------- 1 file changed, 23 insertions(+), 8 deletions(-) diff --git a/meshmode/mesh/generation.py b/meshmode/mesh/generation.py index 65761586..ce2a5f23 100644 --- a/meshmode/mesh/generation.py +++ b/meshmode/mesh/generation.py @@ -175,7 +175,10 @@ def qbx_peanut(t): # {{{ make_curve_mesh -def make_curve_mesh(curve_f, element_boundaries, order): +def make_curve_mesh(curve_f, element_boundaries, order, + unit_nodes=None, + node_vertex_consistency_tolerance=None, + return_parametrization_points=False): """ :arg curve_f: A callable representing a parametrization for a curve, accepting a vector of point locations and returning @@ -183,15 +186,20 @@ def make_curve_mesh(curve_f, element_boundaries, order): :arg element_boundaries: a vector of element boundary locations in :math:`[0,1]`, in order. 0 must be the first entry, 1 the last one. - :returns: a :class:`meshmode.mesh.Mesh` + :arg unit_nodes: if given, the unit nodes to use. Must have shape + ``(dim, nnoodes)``. + :returns: a :class:`meshmode.mesh.Mesh`, or if *return_parametrization_points* + is True, a tuple ``(mesh, par_points)``, where *par_points* is an array of + parametrization points. """ assert element_boundaries[0] == 0 assert element_boundaries[-1] == 1 nelements = len(element_boundaries) - 1 - unodes = mp.warp_and_blend_nodes(1, order) - nodes_01 = 0.5*(unodes+1) + if unit_nodes is None: + unit_nodes = mp.warp_and_blend_nodes(1, order) + nodes_01 = 0.5*(unit_nodes+1) vertices = curve_f(element_boundaries) @@ -200,7 +208,8 @@ def make_curve_mesh(curve_f, element_boundaries, order): # (el_nr, node_nr) t = el_starts[:, np.newaxis] + el_lengths[:, np.newaxis]*nodes_01 - nodes = curve_f(t.ravel()).reshape(vertices.shape[0], nelements, -1) + t = t.ravel() + nodes = curve_f(t).reshape(vertices.shape[0], nelements, -1) from meshmode.mesh import Mesh, SimplexElementGroup egroup = SimplexElementGroup( @@ -210,12 +219,18 @@ def make_curve_mesh(curve_f, element_boundaries, order): np.arange(1, nelements+1, dtype=np.int32) % nelements, ]).T, nodes=nodes, - unit_nodes=unodes) + unit_nodes=unit_nodes) - return Mesh( + mesh = Mesh( vertices=vertices, groups=[egroup], nodal_adjacency=None, - facial_adjacency_groups=None) + facial_adjacency_groups=None, + node_vertex_consistency_tolerance=node_vertex_consistency_tolerance) + + if return_parametrization_points: + return mesh, t + else: + return mesh # }}} -- GitLab From 403082a41ca827e1790eaf4462f7fd94131aebab Mon Sep 17 00:00:00 2001 From: Andreas Kloeckner Date: Fri, 16 Sep 2016 19:03:09 -0500 Subject: [PATCH 19/23] Check length of refine_flags in refiner --- meshmode/mesh/refinement/__init__.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index ac2f9833..cb001e32 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -348,6 +348,10 @@ class Refiner(object): indicating which elements should be split. """ + if len(refine_flags) != self.last_mesh.nelements: + raise ValueError("length of refine_flags does not match " + "element count of last generated mesh") + #vertices and groups for next generation nvertices = len(self.last_mesh.vertices[0]) -- GitLab From 527d3fc9bda539bf7fe5ddf5102601ef47c36f7f Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 19 Sep 2016 22:21:03 -0500 Subject: [PATCH 20/23] Ensure unit_nodes has the right shape for the 1d case of QuadratureSimplexElementGroup --- meshmode/discretization/connection/refinement.py | 2 ++ meshmode/discretization/poly_element.py | 8 +++++++- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/meshmode/discretization/connection/refinement.py b/meshmode/discretization/connection/refinement.py index 4477a9bf..5ef91fe7 100644 --- a/meshmode/discretization/connection/refinement.py +++ b/meshmode/discretization/connection/refinement.py @@ -51,6 +51,8 @@ def _map_unit_nodes_to_children(unit_nodes, tesselation): """ ref_vertices = np.array(tesselation.ref_vertices, dtype=np.float) + assert len(unit_nodes.shape) == 2 + for child_element in tesselation.children: center = np.vstack(ref_vertices[child_element[0]]) # Scale by 1/2 since sides in the tesselation have length 2. diff --git a/meshmode/discretization/poly_element.py b/meshmode/discretization/poly_element.py index a5319bc8..54aaabef 100644 --- a/meshmode/discretization/poly_element.py +++ b/meshmode/discretization/poly_element.py @@ -144,7 +144,13 @@ class QuadratureSimplexElementGroup(PolynomialSimplexElementGroupBase): @property @memoize_method def unit_nodes(self): - return self._quadrature_rule().nodes + result = self._quadrature_rule().nodes + if len(result.shape) == 1: + result = np.array([result]) + + dim2, nunit_nodes = result.shape + assert dim2 == self.mesh_el_group.dim + return result @property @memoize_method -- GitLab From e62f12618dacd567f03eb87042cf95fc84d91adb Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 19 Sep 2016 19:23:25 -0500 Subject: [PATCH 21/23] Refiner: Make 1D curve refinement work. --- meshmode/mesh/refinement/__init__.py | 11 ++++++----- meshmode/mesh/tesselate.py | 6 ++++++ 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/meshmode/mesh/refinement/__init__.py b/meshmode/mesh/refinement/__init__.py index cb001e32..16639442 100644 --- a/meshmode/mesh/refinement/__init__.py +++ b/meshmode/mesh/refinement/__init__.py @@ -73,17 +73,19 @@ class Refiner(object): # {{{ constructor def __init__(self, mesh): - from meshmode.mesh.tesselate import tesselatetet, tesselatetri + from meshmode.mesh.tesselate import tesselateseg, tesselatetet, tesselatetri self.lazy = False self.seen_tuple = {} self.group_refinement_records = [] + seg_node_tuples, seg_result = tesselateseg() tri_node_tuples, tri_result = tesselatetri() tet_node_tuples, tet_result = tesselatetet() #quadrilateral_node_tuples = [ #print tri_result, tet_result - self.simplex_node_tuples = [None, None, tri_node_tuples, tet_node_tuples] + self.simplex_node_tuples = [ + None, seg_node_tuples, tri_node_tuples, tet_node_tuples] # Dimension-parameterized tesselations for refinement - self.simplex_result = [None, None, tri_result, tet_result] + self.simplex_result = [None, seg_result, tri_result, tet_result] #print tri_node_tuples, tri_result #self.simplex_node_tuples, self.simplex_result = tesselatetet() self.last_mesh = mesh @@ -594,8 +596,7 @@ class Refiner(object): midpoint_vertices = [] vertex_indices = grp.vertex_indices[iel_grp] #if simplex - if len(grp.vertex_indices[iel_grp]) == len(self.last_mesh.vertices)+1: - + if len(vertex_indices) == grp.dim + 1: # {{{ Get midpoints for all pairs of vertices for i in range(len(vertex_indices)): diff --git a/meshmode/mesh/tesselate.py b/meshmode/mesh/tesselate.py index 5d0e03ae..c9527802 100644 --- a/meshmode/mesh/tesselate.py +++ b/meshmode/mesh/tesselate.py @@ -6,6 +6,12 @@ def add_tuples(a, b): return tuple(ac+bc for ac, bc in zip(a, b)) +def tesselateseg(): + node_tuples = [(0,), (1,), (2,)] + result = [(0, 1), (1, 2)] + return [node_tuples, result] + + def tesselatetri(): result = [] -- GitLab From d19d6514469ffed23f85e621737aeb4bd1f9e9d5 Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 19 Sep 2016 22:27:00 -0500 Subject: [PATCH 22/23] Add 1D refinement test. --- meshmode/mesh/refinement/utils.py | 2 +- test/test_refinement.py | 11 ++++++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/meshmode/mesh/refinement/utils.py b/meshmode/mesh/refinement/utils.py index 0d0476f3..d34fa14a 100644 --- a/meshmode/mesh/refinement/utils.py +++ b/meshmode/mesh/refinement/utils.py @@ -80,7 +80,7 @@ def check_nodal_adj_against_geometry(mesh, tol=1e-12): nearby_origin_vertex = mesh.vertices[ :, nearby_grp.vertex_indices[nearby_iel][0]] # noqa transformation = np.empty( - (len(mesh.vertices), nvertices_per_element-1)) + (mesh.ambient_dim, mesh.ambient_dim)) vertex_transformed = vertex - nearby_origin_vertex for inearby_vertex_index, nearby_vertex_index in enumerate( diff --git a/test/test_refinement.py b/test/test_refinement.py index 9c36b975..3f395a36 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -31,7 +31,7 @@ from pyopencl.tools import ( # noqa pytest_generate_tests_for_pyopencl as pytest_generate_tests) from meshmode.mesh.generation import ( # noqa - generate_icosahedron, generate_box_mesh) + generate_icosahedron, generate_box_mesh, make_curve_mesh, ellipse) from meshmode.mesh.refinement.utils import check_nodal_adj_against_geometry from meshmode.mesh.refinement import Refiner @@ -80,6 +80,15 @@ def uniform_refine_flags(mesh): # partial(random_refine_flags, 0.4), # 3), + ("3_to_1_ellipse_unif", + partial( + make_curve_mesh, + partial(ellipse, 3), + np.linspace(0, 1, 20), + order=2), + uniform_refine_flags, + 4), + ("rect2d_rand", partial(generate_box_mesh, ( np.linspace(0, 1, 3), -- GitLab From f8166d075f7b946639fb31c27b677b3e639fdead Mon Sep 17 00:00:00 2001 From: Matt Wala Date: Mon, 19 Sep 2016 22:47:35 -0500 Subject: [PATCH 23/23] Another 1d refinement test. --- test/test_refinement.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/test/test_refinement.py b/test/test_refinement.py index 3f395a36..fca3aa93 100644 --- a/test/test_refinement.py +++ b/test/test_refinement.py @@ -147,6 +147,7 @@ def test_refinement(case_name, mesh_gen, flag_gen, num_generations): PolynomialEquidistantGroupFactory ]) @pytest.mark.parametrize(("mesh_name", "dim", "mesh_pars"), [ + ("circle", 1, [20, 30, 40]), ("blob", 2, [1e-1, 8e-2, 5e-2]), ("warp", 2, [4, 5, 6]), ("warp", 3, [4, 5, 6]), @@ -180,7 +181,12 @@ def test_refinement_connection( for mesh_par in mesh_pars: # {{{ get mesh - if mesh_name == "blob": + if mesh_name == "circle": + assert dim == 1 + h = 1 / mesh_par + mesh = make_curve_mesh( + partial(ellipse, 1), np.linspace(0, 1, mesh_par + 1), order=1) + elif mesh_name == "blob": assert dim == 2 h = mesh_par mesh = gen_blob_mesh(h) -- GitLab