Source code for pymor.grids.subgrid

# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)

import weakref

import numpy as np

from pymor.grids.boundaryinfos import GenericBoundaryInfo
from pymor.grids.interfaces import AffineGridInterface


[docs]class SubGrid(AffineGridInterface): """A subgrid of a |Grid|. Given a |Grid| and a list of codim-0 entities we construct the minimal subgrid of the grid, containing all the given entities. Parameters ---------- parent_grid |Grid| of which a subgrid is to be created. parent_entities |NumPy array| of global indices of the codim-0 entities which are to be contained in the subgrid. Attributes ---------- parent_grid The |Grid| from which the subgrid was constructed. :class:`Subgrid` only stores a :mod:`weakref` to the grid, so accessing this property might return `None` if the original grid has been destroyed. """ reference_element = None def __init__(self, parent_grid, parent_entities): assert parent_grid is not None, \ 'parent_grid is None. Maybe you have called sub_grid.with(parent_entities=e)\n' \ 'on a SubGrid for which the parent grid has been destroyed?' assert isinstance(parent_grid, AffineGridInterface) self.dim = parent_grid.dim self.reference_element = parent_grid.reference_element parent_indices = [np.array(np.unique(parent_entities), dtype=np.int32)] assert len(parent_indices[0] == len(parent_entities)) subentities = [np.arange(len(parent_indices[0]), dtype=np.int32).reshape((-1, 1))] for codim in range(1, self.dim + 1): SUBE = parent_grid.subentities(0, codim)[parent_indices[0]] if np.any(SUBE < 0): raise NotImplementedError UI, UI_inv = np.unique(SUBE, return_inverse=True) subentities.append(np.array(UI_inv.reshape(SUBE.shape), dtype=np.int32)) parent_indices.append(np.array(UI, dtype=np.int32)) self.parent_entities = parent_entities self.__parent_grid = weakref.ref(parent_grid) self.__parent_indices = parent_indices self.__subentities = subentities embeddings = parent_grid.embeddings(0) self.__embeddings = (embeddings[0][parent_indices[0]], embeddings[1][parent_indices[0]]) @property def parent_grid(self): return self.__parent_grid()
[docs] def parent_indices(self, codim): """`retval[e]` is the index of the `e`-th codim-`codim` entity in the parent grid.""" assert 0 <= codim <= self.dim, 'Invalid codimension' return self.__parent_indices[codim]
[docs] def indices_from_parent_indices(self, ind, codim): """Maps a |NumPy array| of indicies of codim-`codim` entites of the parent grid to indicies of the subgrid. Raises ------ ValueError Not all provided indices correspond to entities contained in the subgrid. """ assert 0 <= codim <= self.dim, 'Invalid codimension' ind = ind.ravel() # TODO Find better implementation of the following R = np.argmax(ind[:, np.newaxis] - self.__parent_indices[codim][np.newaxis, :] == 0, axis=1) if not np.all(self.__parent_indices[codim][R] == ind): raise ValueError('Not all parent indices found') return np.array(R, dtype=np.int32)
[docs] def size(self, codim): assert 0 <= codim <= self.dim, 'Invalid codimension' return len(self.__parent_indices[codim])
[docs] def subentities(self, codim, subentity_codim): if codim == 0: assert codim <= subentity_codim <= self.dim, 'Invalid subentity codimension' return self.__subentities[subentity_codim] else: return super().subentities(codim, subentity_codim)
[docs] def embeddings(self, codim): if codim == 0: return self.__embeddings else: return super().embeddings(codim)
def __getstate__(self): d = self.__dict__.copy() del d['_SubGrid__parent_grid'] return d
[docs]def make_sub_grid_boundary_info(sub_grid, parent_grid, parent_grid_boundary_info, new_boundary_type=None): """Derives a |BoundaryInfo| for a :class:`~pymor.grids.subgrid.SubGrid`. Parameters ---------- sub_grid The :class:`~pymor.grids.subgrid.SubGrid` for which a |BoundaryInfo| is created. parent_grid The parent |Grid|. parent_grid_boundary_info The |BoundaryInfo| of the parent |Grid| from which to derive the |BoundaryInfo| new_boundary_type The boundary type which is assigned to the new boundaries of `subgrid`. If `None`, no boundary type is assigned. Returns ------- |BoundaryInfo| associated with sub_grid. """ boundary_types = parent_grid_boundary_info.boundary_types masks = {} for codim in range(1, sub_grid.dim + 1): parent_indices = sub_grid.parent_indices(codim)[sub_grid.boundaries(codim)] new_boundaries = np.where(np.logical_not(parent_grid.boundary_mask(codim)[parent_indices])) new_boundaries_sg_indices = sub_grid.boundaries(codim)[new_boundaries] for t in boundary_types: m = parent_grid_boundary_info.mask(t, codim)[sub_grid.parent_indices(codim)] if t == new_boundary_type: m[new_boundaries_sg_indices] = True if codim == 1: masks[t] = [m] else: masks[t].append(m) if new_boundary_type is not None and new_boundary_type not in boundary_types: m = np.zeros(sub_grid.size(codim), dtype=np.bool) m[new_boundaries_sg_indices] = True if codim == 1: masks[new_boundary_type] = [m] else: masks[new_boundary_type].append(m) return GenericBoundaryInfo(sub_grid, masks)