# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2020 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
import numpy as np
from pymor.discretizers.builtin.grids.interfaces import AffineGridWithOrthogonalCenters
from pymor.discretizers.builtin.grids.referenceelements import square
[docs]class RectGrid(AffineGridWithOrthogonalCenters):
"""Basic implementation of a rectangular |Grid| on a rectangular domain.
The global face, edge and vertex indices are given as follows ::
x1
^
|
6--10---7--11---8
| | |
3 2 4 3 5
| | |
3---8---4---9---5
| | |
0 0 1 1 2
| | |
0---6---1---7---2 --> x0
Parameters
----------
num_intervals
Tuple `(n0, n1)` determining a grid with `n0` x `n1` codim-0 entities.
domain
Tuple `(ll, ur)` where `ll` defines the lower left and `ur` the upper right
corner of the domain.
identify_left_right
If `True`, the left and right boundaries are identified, i.e. the left-most
codim-0 entities become neighbors of the right-most codim-0 entities.
identify_bottom_top
If `True`, the bottom and top boundaries are identified, i.e. the bottom-most
codim-0 entities become neighbors of the top-most codim-0 entities.
"""
dim = 2
reference_element = square
def __init__(self, num_intervals=(2, 2), domain=([0, 0], [1, 1]),
identify_left_right=False, identify_bottom_top=False):
if identify_left_right:
assert num_intervals[0] > 1
if identify_bottom_top:
assert num_intervals[1] > 1
domain = np.array(domain)
self.__auto_init(locals())
self.x0_num_intervals = num_intervals[0]
self.x1_num_intervals = num_intervals[1]
self.x0_range = self.domain[:, 0]
self.x1_range = self.domain[:, 1]
self.x0_width = self.x0_range[1] - self.x0_range[0]
self.x1_width = self.x1_range[1] - self.x1_range[0]
self.x0_diameter = self.x0_width / self.x0_num_intervals
self.x1_diameter = self.x1_width / self.x1_num_intervals
self.diameter_max = max(self.x0_diameter, self.x1_diameter)
self.diameter_min = min(self.x0_diameter, self.x1_diameter)
n_elements = self.x0_num_intervals * self.x1_num_intervals
ni0, ni1 = num_intervals
# TOPOLOGY
# mapping of structured indices to global indices
structured_to_global_0 = np.arange(ni0 * ni1, dtype=np.int32).reshape((ni1, ni0)).swapaxes(0, 1)
structured_to_global_2 = (np.arange((ni0 + 1) * (ni1 + 1), dtype=np.int32)
.reshape(((ni1 + 1), (ni0 + 1))).swapaxes(0, 1))
self._structured_to_global = [structured_to_global_0, None, structured_to_global_2]
global_to_structured_0 = np.empty((ni0 * ni1, 2), dtype=np.int32)
global_to_structured_0[:, 0] = np.tile(np.arange(ni0, dtype=np.int32), ni1)
global_to_structured_0[:, 1] = np.repeat(np.arange(ni1, dtype=np.int32), ni0)
global_to_structured_1 = np.empty(((ni0 + 1) * (ni1 + 1), 2), dtype=np.int32)
global_to_structured_1[:, 0] = np.tile(np.arange(ni0 + 1, dtype=np.int32), ni1 + 1)
global_to_structured_1[:, 1] = np.repeat(np.arange(ni1 + 1, dtype=np.int32), ni0 + 1)
self._global_to_structured = [global_to_structured_0, None, global_to_structured_1]
# calculate subentities -- codim-0
codim1_subentities = np.empty((ni1, ni0, 4), dtype=np.int32)
if identify_left_right:
codim1_subentities[:, :, 3] = np.arange(ni0 * ni1).reshape((ni1, ni0))
codim1_subentities[:, :, 1] = codim1_subentities[:, :, 3] + 1
codim1_subentities[:, -1, 1] = codim1_subentities[:, 0, 3]
else:
codim1_subentities[:, :, 3] = (np.arange(ni0)[np.newaxis, :]
+ (np.arange(ni1) * (ni0 + 1))[:, np.newaxis])
codim1_subentities[:, :, 1] = codim1_subentities[:, :, 3] + 1
offset = np.max(codim1_subentities[:, :, [1, 3]]) + 1
codim1_subentities[:, :, 0] = (np.arange(ni0 * ni1) + offset).reshape((ni1, ni0))
codim1_subentities[:, :, 2] = codim1_subentities[:, :, 0] + num_intervals[0]
if identify_bottom_top:
codim1_subentities[-1, :, 2] = codim1_subentities[0, :, 0]
codim1_subentities = codim1_subentities.reshape((-1, 4))
# calculate subentities -- codim-1
codim2_subentities = np.empty((ni1, ni0, 4), dtype=np.int32)
if identify_left_right:
codim2_subentities[:, :, 0] = np.arange(ni0 * ni1).reshape((ni1, ni0))
codim2_subentities[:, :, 1] = codim2_subentities[:, :, 0] + 1
codim2_subentities[:, -1, 1] = codim2_subentities[:, 0, 0]
codim2_subentities[:, :, 3] = codim2_subentities[:, :, 0] + ni0
codim2_subentities[:, :, 2] = codim2_subentities[:, :, 1] + ni0
else:
codim2_subentities[:, :, 0] = (np.arange(ni0)[np.newaxis, :]
+ (np.arange(ni1) * (ni0 + 1))[:, np.newaxis])
codim2_subentities[:, :, 1] = codim2_subentities[:, :, 0] + 1
codim2_subentities[:, :, 3] = codim2_subentities[:, :, 0] + (ni0 + 1)
codim2_subentities[:, :, 2] = codim2_subentities[:, :, 0] + (ni0 + 2)
if identify_bottom_top:
codim2_subentities[-1, :, 3] = codim2_subentities[0, :, 0]
codim2_subentities[-1, :, 2] = codim2_subentities[0, :, 1]
codim2_subentities = codim2_subentities.reshape((-1, 4))
self.__subentities = (codim1_subentities, codim2_subentities)
# calculate number of elements
self.__sizes = (n_elements,
np.max(codim1_subentities) + 1,
np.max(codim2_subentities) + 1)
# GEOMETRY
# embeddings
x0_shifts = np.arange(self.x0_num_intervals) * self.x0_diameter + self.x0_range[0]
x1_shifts = np.arange(self.x1_num_intervals) * self.x1_diameter + self.x1_range[0]
shifts = np.array(np.meshgrid(x0_shifts, x1_shifts)).reshape((2, -1))
A = np.tile(np.diag([self.x0_diameter, self.x1_diameter]), (n_elements, 1, 1))
B = shifts.T
self.__embeddings = (A, B)
[docs] def __reduce__(self):
return (RectGrid,
(self.num_intervals, self.domain, self.identify_left_right, self.identify_bottom_top))
[docs] def __str__(self):
return (f'Rect-Grid on domain '
f'[{self.x0_range[0]},{self.x0_range[1]}] x [{self.x1_range[0]},{self.x1_range[1]}]\n'
f'x0-intervals: {self.x0_num_intervals}, x1-intervals: {self.x1_num_intervals}\n'
f'faces: {self.size(0)}, edges: {self.size(1)}, vertices: {self.size(2)}')
[docs] def size(self, codim=0):
assert 0 <= codim <= 2, 'Invalid codimension'
return self.__sizes[codim]
[docs] def subentities(self, codim, subentity_codim):
assert 0 <= codim <= 2, 'Invalid codimension'
assert codim <= subentity_codim <= self.dim, 'Invalid subentity codimension'
if codim == 0:
if subentity_codim == 0:
return np.arange(self.size(0), dtype='int32')[:, np.newaxis]
else:
return self.__subentities[subentity_codim - 1]
else:
return super().subentities(codim, subentity_codim)
[docs] def embeddings(self, codim=0):
if codim == 0:
return self.__embeddings
else:
return super().embeddings(codim)
[docs] def bounding_box(self):
return np.array(self.domain)
[docs] def structured_to_global(self, codim):
"""Returns an |array| which maps structured indices to global codim-`codim` indices.
In other words `structured_to_global(codim)[i, j]` is the global index of the i-th in
x0-direction and j-th in x1-direction codim-`codim` entity of the grid.
"""
if self.identify_left_right or self.identify_bottom_top or codim not in (0, 2):
raise NotImplementedError
return self._structured_to_global[codim]
[docs] def global_to_structured(self, codim):
"""Returns an array which maps global codim-`codim` indices to structured indices.
I.e. if `GTS = global_to_structured(codim)` and `STG = structured_to_global(codim)`, then
`STG[GTS[:, 0], GTS[:, 1]] == numpy.arange(size(codim))`.
"""
if self.identify_left_right or self.identify_bottom_top or codim not in (0, 2):
raise NotImplementedError
return self._global_to_structured[codim]
[docs] def vertex_coordinates(self, dim):
"""Returns an array of the x_dim coordinates of the grid vertices.
I.e. ::
centers(2)[structured_to_global(2)[i, j]] == np.array([vertex_coordinates(0)[i], vertex_coordinates(1)[j]])
"""
assert 0 <= dim < 2
return np.linspace(self.domain[0, dim], self.domain[1, dim], self.num_intervals[dim] + 1)
[docs] def orthogonal_centers(self):
return self.centers(0)
[docs] def visualize(self, U, codim=2, **kwargs):
"""Visualize scalar data associated to the grid as a patch plot.
Parameters
----------
U
|NumPy array| of the data to visualize. If `U.dim == 2 and len(U) > 1`, the
data is visualized as a time series of plots. Alternatively, a tuple of
|Numpy arrays| 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.discretizers.builtin.gui.qt.visualize_patch`
"""
from pymor.discretizers.builtin.gui.qt import visualize_patch
from pymor.vectorarrays.interface import VectorArray
from pymor.vectorarrays.numpy import NumpyVectorSpace, NumpyVectorArray
if isinstance(U, (np.ndarray, VectorArray)):
U = (U,)
assert all(isinstance(u, (np.ndarray, VectorArray)) for u in U)
U = tuple(NumpyVectorSpace.make_array(u) if isinstance(u, np.ndarray) else
u if isinstance(u, NumpyVectorArray) else
NumpyVectorSpace.make_array(u.to_numpy())
for u in U)
bounding_box = kwargs.pop('bounding_box', self.domain)
visualize_patch(self, U, codim=codim, bounding_box=bounding_box, **kwargs)