Add UnstructuredTriangleGrid

setup.py

@@ -126,7 +126,8 @@ def _setup(**kwargs):
                 'test': PyTest}
     from numpy import get_include
     ext_modules = [Extension("pymor.tools.relations", ["src/pymor/tools/relations.pyx"], include_dirs=[get_include()]),
-                   Extension("pymor.tools.inplace", ["src/pymor/tools/inplace.pyx"], include_dirs=[get_include()])]
+                   Extension("pymor.tools.inplace", ["src/pymor/tools/inplace.pyx"], include_dirs=[get_include()]),
+                   Extension("pymor.grids._unstructured", ["src/pymor/grids/_unstructured.pyx"], include_dirs=[get_include()])]
     # for some reason the *pyx files don't end up in sdist tarballs -> manually add them as package data
     # this _still_ doesn't make them end up in the tarball however -> manually add them in MANIFEST.in
     kwargs['package_data'] = {'pymor': list(itertools.chain(*[i.sources for i in ext_modules])) }

src/pymor/discretizers/advection.py

@@ -119,7 +119,7 @@ def initial_projection(U, mu):
 
     products = {'l2': L2Product(grid, boundary_info)}
     if grid.dim == 2:
-        visualizer = PatchVisualizer(grid=grid, bounding_box=grid.domain, codim=0)
+        visualizer = PatchVisualizer(grid=grid, bounding_box=grid.bounding_box(), codim=0)
     elif grid.dim == 1:
         visualizer = Matplotlib1DVisualizer(grid, codim=0)
     else:

src/pymor/discretizers/elliptic.py

@@ -11,9 +11,7 @@
 from pymor.discretizations.basic import StationaryDiscretization
 from pymor.domaindiscretizers.default import discretize_domain_default
 from pymor.grids.boundaryinfos import EmptyBoundaryInfo
-from pymor.grids.oned import OnedGrid
-from pymor.grids.rect import RectGrid
-from pymor.grids.tria import TriaGrid
+from pymor.grids.referenceelements import line, triangle, square
 from pymor.gui.qt import PatchVisualizer, Matplotlib1DVisualizer
 from pymor.operators import cg, fv
 from pymor.operators.constructions import LincombOperator
@@ -64,9 +62,9 @@ def discretize_elliptic_cg(analytical_problem, diameter=None, domain_discretizer
         else:
             grid, boundary_info = domain_discretizer(analytical_problem.domain, diameter=diameter)
 
-    assert isinstance(grid, (OnedGrid, TriaGrid, RectGrid))
+    assert grid.reference_element in (line, triangle, square)
 
-    if isinstance(grid, RectGrid):
+    if grid.reference_element is square:
         Operator = cg.DiffusionOperatorQ1
         Functional = cg.L2ProductFunctionalQ1
     else:
@@ -91,13 +89,15 @@ def discretize_elliptic_cg(analytical_problem, diameter=None, domain_discretizer
 
     F = Functional(grid, p.rhs, boundary_info, dirichlet_data=p.dirichlet_data, neumann_data=p.neumann_data)
 
-    if isinstance(grid, (TriaGrid, RectGrid)):
-        visualizer = PatchVisualizer(grid=grid, bounding_box=grid.domain, codim=2)
-    else:
+    if grid.reference_element in (triangle, square):
+        visualizer = PatchVisualizer(grid=grid, bounding_box=grid.bounding_box(), codim=2)
+    elif grid.reference_element is line:
         visualizer = Matplotlib1DVisualizer(grid=grid, codim=1)
+    else:
+        visualizer = None
 
     empty_bi = EmptyBoundaryInfo(grid)
-    l2_product = cg.L2ProductQ1(grid, empty_bi) if isinstance(grid, RectGrid) else cg.L2ProductP1(grid, empty_bi)
+    l2_product = cg.L2ProductQ1(grid, empty_bi) if grid.reference_element is square else cg.L2ProductP1(grid, empty_bi)
     h1_semi_product = Operator(grid, empty_bi)
     products = {'h1': l2_product + h1_semi_product,
                 'h1_semi': h1_semi_product,
@@ -182,9 +182,9 @@ def discretize_elliptic_fv(analytical_problem, diameter=None, domain_discretizer
         F = fv.L2ProductFunctional(grid, p.rhs, boundary_info=boundary_info, dirichlet_data=p.dirichlet_data,
                                    diffusion_function=p.diffusion_functions[0], neumann_data=p.neumann_data)
 
-    if isinstance(grid, (TriaGrid, RectGrid)):
-        visualizer = PatchVisualizer(grid=grid, bounding_box=grid.domain, codim=0)
-    elif isinstance(grid, (OnedGrid)):
+    if grid.reference_element in (triangle, square):
+        visualizer = PatchVisualizer(grid=grid, bounding_box=grid.bounding_box(), codim=0)
+    elif grid.reference_element is line:
         visualizer = Matplotlib1DVisualizer(grid=grid, codim=0)
     else:
         visualizer = None

src/pymor/grids/_unstructured.pyx

@@ -0,0 +1,74 @@
+# This file is part of the pyMOR project (http://www.pymor.org).
+# Copyright Holders: Rene Milk, Stephan Rave, Felix Schindler
+# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+
+from __future__ import division
+
+import numpy as np
+import cython
+cimport numpy as np
+
+DTYPE = np.int32
+ctypedef np.int32_t DTYPE_t
+
+
+@cython.boundscheck(False)
+def compute_edges(np.ndarray[DTYPE_t, ndim=2] faces, int num_vertices):
+    cdef unsigned int i
+    cdef unsigned int j
+    cdef unsigned int m
+    cdef unsigned int n
+    cdef int x
+    cdef unsigned int index
+    cdef unsigned int max_edges
+
+    cdef np.ndarray[DTYPE_t, ndim=1] edge_counts
+    edge_counts = np.zeros(num_vertices, dtype=DTYPE)
+
+    for i in xrange(faces.shape[0]):
+        for j in xrange(faces.shape[1]):
+            edge_counts[<unsigned int>faces[i, j]] += 2
+
+    max_edges = np.max(edge_counts)
+    del edge_counts
+
+    cdef unsigned int edge_counter
+    cdef np.ndarray[DTYPE_t, ndim=2] edge_memo_vertices
+    cdef np.ndarray[DTYPE_t, ndim=2] edge_memo_ids
+    cdef np.ndarray[DTYPE_t, ndim=2] edges
+    edge_memo_vertices = np.empty((num_vertices, max_edges), dtype=DTYPE)
+    edge_memo_vertices[:] = -1
+    edge_memo_ids= np.empty((num_vertices, max_edges), dtype=DTYPE)
+    edges = np.empty_like(faces)
+
+    edge_counter = 0
+    for i in xrange(faces.shape[0]):
+        for j in xrange(3):
+            if j == 0:
+                m = faces[i, <unsigned int>1]
+                n = faces[i, <unsigned int>2]
+            elif j == 1:
+                m = faces[i, <unsigned int>2]
+                n = faces[i, <unsigned int>0]
+            else:
+                m = faces[i, <unsigned int>0]
+                n = faces[i, <unsigned int>1]
+
+            if m > n:
+                m, n = n, m
+
+            # try to find edge in memo
+            for index in xrange(max_edges):
+                x = edge_memo_vertices[m, index]
+                if x == n or x == -1:
+                    break
+
+            if x == n:
+                edges[i, j] = edge_memo_ids[m, index]
+            else:
+                assert x == -1
+                edge_memo_ids[m, index] = edges[i, j] = edge_counter
+                edge_memo_vertices[m, index] = n
+                edge_counter += 1
+
+    return edges, edge_counter

src/pymor/grids/defaultimpl.py

@@ -265,3 +265,12 @@ def _quadrature_points(self, codim, order, npoints, quadrature_type):
         P, _ = self.reference_element(codim).quadrature(order, npoints, quadrature_type)
         A, B = self.embeddings(codim)
         return np.einsum('eij,kj->eki', A, P) + B[:, np.newaxis, :]
+
+    @cached
+    def _bounding_box(self):
+        bbox = np.empty((2, self.dim_outer))
+        centers = self.centers(self.dim)
+        for dim in range(self.dim_outer):
+            bbox[0, dim] = np.min(centers[:, dim])
+            bbox[1, dim] = np.max(centers[:, dim])
+        return bbox

src/pymor/grids/interfaces.py

@@ -291,6 +291,10 @@ def quadrature_points(self, codim, order=None, npoints=None, quadrature_type='de
         """
         return self._quadrature_points(codim, order, npoints, quadrature_type)
 
+    def bounding_box(self):
+        """returns a `(2, dim_outer)`-shaped array containing lower/upper bounding box coordinates."""
+        return self._bounding_box()
+
 
 class AffineGridWithOrthogonalCentersInterface(AffineGridInterface):
     """|AffineGrid| with an additional `orthogonal_centers` method."""

src/pymor/grids/oned.py

@@ -27,6 +27,7 @@ class OnedGrid(AffineGridWithOrthogonalCentersInterface):
     reference_element = line
 
     def __init__(self, domain=(0, 1), num_intervals=4, identify_left_right=False):
+        assert domain[0] < domain[1]
         self.reference_element = line
         self._domain = np.array(domain)
         self._num_intervals = num_intervals
@@ -74,6 +75,9 @@ def embeddings(self, codim):
         else:
             return super(OnedGrid, self).embeddings(codim)
 
+    def bounding_box(self):
+        return np.array(self._domain).reshape((2, 1))
+
     def orthogonal_centers(self):
         return self.centers(0)
 

src/pymor/grids/rect.py

@@ -167,6 +167,9 @@ def embeddings(self, codim=0):
         else:
             return super(RectGrid, self).embeddings(codim)
 
+    def bounding_box(self):
+        return np.array(self.domain)
+
     def structured_to_global(self, codim):
         """Returns an array which maps structured indices to global codim-`codim` indices.
 

src/pymor/grids/tria.py

@@ -197,6 +197,9 @@ def embeddings(self, codim=0):
         else:
             return super(TriaGrid, self).embeddings(codim)
 
+    def bounding_box(self):
+        return np.array(self.domain)
+
     @cached
     def orthogonal_centers(self):
         embeddings = self.embeddings(0)

src/pymor/grids/unstructured.py

@@ -0,0 +1,95 @@
+# This file is part of the pyMOR project (http://www.pymor.org).
+# Copyright Holders: Rene Milk, Stephan Rave, Felix Schindler
+# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+
+from pymor.grids.interfaces import AffineGridInterface
+from pymor.grids.referenceelements import triangle
+from pymor.grids._unstructured import compute_edges
+
+
+class UnstructuredTriangleGrid(AffineGridInterface):
+    """A generic unstructured, triangular grid.
+
+    Parameters
+    ----------
+    vertices
+        A (num_vertices, 2)-shaped |array| containing the coordinates
+        of all vertices in the grid. The row numbers in the array will
+        be the global indices of the given vertices (codim 2 entities).
+    faces
+        A (num_faces, 3)-shaped |array| containing the global indices
+        of the vertices which define a given triangle in the grid.
+        The row numbers in the array will be the global indices of the
+        given triangles (codim 0 entities).
+    """
+
+    dim = 2
+    dim_outer = 2
+    reference_element = triangle
+
+    def __init__(self, vertices, faces):
+        assert faces.shape[1] == 3
+        assert np.min(faces) == 0
+        assert np.max(faces) == len(vertices) - 1
+
+        vertices = vertices.astype(np.float64, copy=False)
+        faces = faces.astype(np.int32, copy=False)
+        edges, num_edges = compute_edges(faces, len(vertices))
+
+        COORDS = vertices[faces]
+        SHIFTS = COORDS[:, 0, :]
+        TRANS = COORDS[:, 1:, :] - SHIFTS[:, np.newaxis, :]
+        TRANS = TRANS.swapaxes(1, 2)
+
+        self.__embeddings = (TRANS, SHIFTS)
+        self.__subentities = (np.arange(len(faces), dtype=np.int32).reshape(-1, 1), edges, faces)
+        self.__sizes = (len(faces), num_edges, len(vertices))
+
+    def size(self, codim=0):
+        assert 0 <= codim <= 2, 'Invalid codimension'
+        return self.__sizes[codim]
+
+    def subentities(self, codim=0, subentity_codim=None):
+        assert 0 <= codim <= 2, 'Invalid codimension'
+        if subentity_codim is None:
+            subentity_codim = codim + 1
+        assert codim <= subentity_codim <= self.dim, 'Invalid subentity codimensoin'
+        if codim == 0:
+            return self.__subentities[subentity_codim]
+        else:
+            return super(UnstructuredTriangleGrid, self).subentities(codim, subentity_codim)
+
+    def embeddings(self, codim=0):
+        if codim == 0:
+            return self.__embeddings
+        else:
+            return super(UnstructuredTriangleGrid, self).embeddings(codim)
+
+    def visualize(self, U, codim=2, **kwargs):
+        """Visualize scalar data associated to the grid as a patch plot.
+
+        Parameters
+        ----------
+        U
+            |VectorArray| of the data to visualize. If `len(U) > 1`, the data is visualized
+            as a time series of plots. Alternatively, a tuple of |VectorArrays| can be
+            provided, in which case a subplot is created for each entry of the tuple. The
+            lengths of all arrays have to agree.
+        codim
+            The codimension of the entities the data in `U` is attached to (either 0 or 2).
+        kwargs
+            See :func:`~pymor.gui.qt.visualize_patch`
+        """
+        from pymor.gui.qt import visualize_patch
+        from pymor.vectorarrays.numpy import NumpyVectorArray
+        if not isinstance(U, NumpyVectorArray):
+            U = NumpyVectorArray(U, copy=False)
+        bounding_box = kwargs.pop('bounding_box', self.bounding_box())
+        visualize_patch(self, U, codim=codim, bounding_box=bounding_box, **kwargs)
+
+    def __str__(self):
+        return 'UnstructuredTriangleGrid with {} triangles, {} edges, {} vertices'.format(*self.__sizes)

src/pymor/gui/qt.py

@@ -34,8 +34,7 @@
 from pymor.core.interfaces import BasicInterface
 from pymor.core.logger import getLogger
 from pymor.grids.oned import OnedGrid
-from pymor.grids.rect import RectGrid
-from pymor.grids.tria import TriaGrid
+from pymor.grids.referenceelements import triangle, square
 from pymor.gui.gl import GLPatchWidget, ColorBarWidget, HAVE_GL
 from pymor.gui.matplotlib import Matplotlib1DWidget, MatplotlibPatchWidget, HAVE_MATPLOTLIB
 from pymor.tools.vtkio import HAVE_PYVTK, write_vtk
@@ -465,7 +464,8 @@ class PatchVisualizer(BasicInterface):
     """
 
     def __init__(self, grid, bounding_box=([0, 0], [1, 1]), codim=2, backend=None, block=False):
-        assert isinstance(grid, (RectGrid, TriaGrid))
+        assert grid.reference_element in (triangle, square)
+        assert grid.dim_outer == 2
         assert codim in (0, 2)
         self.grid = grid
         self.bounding_box = bounding_box

src/pymortests/affine_grid.py

@@ -294,6 +294,15 @@ def test_quadrature_points_values(grid):
                 np.testing.assert_allclose(Q, B[:, np.newaxis, :] + np.einsum('eij,qj->eqi', A, q))
 
 
+def test_bounding_box(grid):
+    g = grid
+    bbox = g.bounding_box()
+    assert bbox.shape == (2, g.dim_outer)
+    assert np.all(bbox[0] <= bbox[1])
+    assert np.all(g.centers(g.dim) >= bbox[0])
+    assert np.all(g.centers(g.dim) <= bbox[1])
+
+
 def test_orthogonal_centers(grid_with_orthogonal_centers):
     g = grid_with_orthogonal_centers
     C = g.orthogonal_centers()

src/pymortests/fixtures/grid.py

@@ -14,6 +14,7 @@
 from pymor.grids.rect import RectGrid
 from pymor.grids.subgrid import SubGrid
 from pymor.grids.tria import TriaGrid
+from pymor.grids.unstructured import UnstructuredTriangleGrid
 
 
 rect_grid_generators = [lambda arg=arg, kwargs=kwargs: RectGrid(arg, **kwargs) for arg, kwargs in
@@ -48,12 +49,15 @@
 
 oned_grid_generators = [lambda kwargs=kwargs: OnedGrid(**kwargs) for kwargs in
                         [dict(domain=np.array((-2, 2)), num_intervals=10),
-                         dict(domain=np.array((-2, -4)), num_intervals=100),
-                         dict(domain=np.array((-2, -4)), num_intervals=100, identify_left_right=True),
-                         dict(domain=np.array((3, 2)), num_intervals=10),
-                         dict(domain=np.array((3, 2)), num_intervals=10, identify_left_right=True),
+                         dict(domain=np.array((-4, -2)), num_intervals=100),
+                         dict(domain=np.array((-4, -2)), num_intervals=100, identify_left_right=True),
+                         dict(domain=np.array((2, 3)), num_intervals=10),
+                         dict(domain=np.array((2, 3)), num_intervals=10, identify_left_right=True),
                          dict(domain=np.array((1, 2)), num_intervals=10000)]]
 
+unstructured_grid_generators = [lambda: UnstructuredTriangleGrid(np.array([[0, 0], [-1, -1], [1, -1], [1, 1], [-1, 1]]),
+                                                                 np.array([[0, 1, 2], [0, 3, 4], [0, 4, 1]]))]
+
 
 def subgrid_factory(grid_generator, neq, seed):
     np.random.seed(seed)
@@ -79,7 +83,8 @@ def subgrid_factory(grid_generator, neq, seed):
                        (lambda: TriaGrid((24, 24)), 4, 123)]]
 
 
-@pytest.fixture(params=(rect_grid_generators + tria_grid_generators + oned_grid_generators + subgrid_generators))
+@pytest.fixture(params=(rect_grid_generators + tria_grid_generators + oned_grid_generators + subgrid_generators
+                        + unstructured_grid_generators))
 def grid(request):
     return request.param()