Source code for pymor.grids.defaultimpl

# 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 numpy as np

from pymor.core.cache import cached
from pymor.tools.inverse import inv_transposed_two_by_two
from pymor.tools.relations import inverse_relation


[docs]class ConformalTopologicalGridDefaultImplementations: """Provides default informations for |ConformalTopologicalGrids|.""" @cached def _subentities(self, codim, subentity_codim): assert 0 <= codim < self.dim, 'Invalid codimension' if subentity_codim > codim + 1: SE = self.subentities(codim, subentity_codim - 1) SESE = self.subentities(subentity_codim - 1, subentity_codim) # we assume that there is only one geometry type ... num_subsubentities = np.unique(SESE[SE[0]]).size SSE = np.empty((SE.shape[0], num_subsubentities), dtype=np.int32) SSE.fill(-1) for ei in range(SE.shape[0]): X = SESE[SE[ei]].ravel() SSE[ei] = X[np.sort(np.unique(X, return_index=True)[1])] return SSE else: raise NotImplementedError @cached def _superentities_with_indices(self, codim, superentity_codim): assert 0 <= codim <= self.dim, f'Invalid codimension (was {codim})' assert 0 <= superentity_codim <= codim, f'Invalid codimension (was {superentity_codim})' SE = self.subentities(superentity_codim, codim) return inverse_relation(SE, size_rhs=self.size(codim), with_indices=True) @cached def _superentities(self, codim, superentity_codim): return self._superentities_with_indices(codim, superentity_codim)[0] @cached def _superentity_indices(self, codim, superentity_codim): return self._superentities_with_indices(codim, superentity_codim)[1] @cached def _neighbours(self, codim, neighbour_codim, intersection_codim): assert 0 <= codim <= self.dim, 'Invalid codimension' assert 0 <= neighbour_codim <= self.dim, 'Invalid codimension' if intersection_codim is None: if codim == neighbour_codim: intersection_codim = codim + 1 else: intersection_codim = min(codim, neighbour_codim) assert max(codim, neighbour_codim) <= intersection_codim <= self.dim, 'Invalid codimension' if intersection_codim == min(codim, neighbour_codim): if codim < neighbour_codim: return self.subentities(codim, neighbour_codim) elif codim > neighbour_codim: return self.superentities(codim, neighbour_codim) else: return np.zeros((self.size(codim), 0), dtype=np.int32) else: EI = self.subentities(codim, intersection_codim) ISE = self.superentities(intersection_codim, neighbour_codim) NB = np.empty((EI.shape[0], EI.shape[1] * ISE.shape[1]), dtype=np.int32) NB.fill(-1) NB_COUNTS = np.zeros(EI.shape[0], dtype=np.int32) if codim == neighbour_codim: for ii, i in np.ndenumerate(EI): if i >= 0: for _, n in np.ndenumerate(ISE[i]): if n != ii[0] and n not in NB[ii[0]]: NB[ii[0], NB_COUNTS[ii[0]]] = n NB_COUNTS[ii[0]] += 1 else: for ii, i in np.ndenumerate(EI): if i >= 0: for _, n in np.ndenumerate(ISE[i]): if n not in NB[ii[0]]: NB[ii[0], NB_COUNTS[ii[0]]] = n NB_COUNTS[ii[0]] += 1 NB = NB[:NB.shape[0], :NB_COUNTS.max()] return NB @cached def _boundaries(self, codim): assert 0 <= codim <= self.dim, 'Invalid codimension' if codim == 1: SE = self.superentities(1, 0) # a codim-1 entity can have at most 2 superentities, and it is a boundary # if it has only one superentity if SE.shape[1] > 1: return np.where(np.any(SE == -1, axis=1))[0].astype('int32') else: return np.arange(SE.shape[0], dtype='int32') elif codim == 0: B1 = self.boundaries(1) if B1.size > 0: B0 = np.unique(self.superentities(1, 0)[B1]) return B0[1:] if B0[0] == -1 else B0 else: return np.array([], dtype=np.int32) else: B1 = self.boundaries(1) if B1.size > 0: BC = np.unique(self.subentities(1, codim)[B1]) return BC[1:] if BC[0] == -1 else BC else: return np.array([], dtype=np.int32) @cached def _boundary_mask(self, codim): M = np.zeros(self.size(codim), dtype='bool') B = self.boundaries(codim) if B.size > 0: M[self.boundaries(codim)] = True return M
[docs]class ReferenceElementDefaultImplementations: """Provides default implementations for |ReferenceElements|.""" @cached def _subentity_embedding(self, subentity_codim): if subentity_codim > 1: A = [] B = [] for i in range(self.size(subentity_codim)): P = np.where(self.subentities(subentity_codim - 1, subentity_codim) == i) parent_index, local_index = P[0][0], P[1][0] A0, B0 = self.subentity_embedding(subentity_codim - 1) A0 = A0[parent_index] B0 = B0[parent_index] A1, B1 = self.sub_reference_element(subentity_codim - 1).subentity_embedding(1) A1 = A1[local_index] B1 = B1[local_index] A.append(np.dot(A0, A1)) B.append(np.dot(A0, B1) + B0) return np.array(A), np.array(B) else: raise NotImplementedError @cached def _sub_reference_element(self, codim): if codim > 1: return self.sub_reference_element(1).sub_reference_element(codim - 1) else: raise NotImplementedError
[docs]class AffineGridDefaultImplementations: """Provides default implementations for |AffineGrids|.""" @cached def _subentities(self, codim, subentity_codim): assert 0 <= codim <= self.dim, 'Invalid codimension' assert 0 < codim, 'Not implemented' P = self.superentities(codim, codim - 1)[:, 0] # we assume here that superentities() is sorted by global index I = self.superentity_indices(codim, codim - 1)[:, 0] SE = self.subentities(codim - 1, subentity_codim)[P] RSE = self.reference_element(codim - 1).subentities(1, subentity_codim - (codim - 1))[I] SSE = np.empty_like(RSE) for i in range(RSE.shape[0]): SSE[i, :] = SE[i, RSE[i]] return SSE @cached def _embeddings(self, codim): assert codim > 0, NotImplemented E = self.superentities(codim, codim - 1)[:, 0] I = self.superentity_indices(codim, codim - 1)[:, 0] A0, B0 = self.embeddings(codim - 1) A0 = A0[E] B0 = B0[E] A1, B1 = self.reference_element(codim - 1).subentity_embedding(1) A = np.zeros((E.shape[0], A0.shape[1], A1.shape[2])) B = np.zeros((E.shape[0], A0.shape[1])) for i in range(A1.shape[0]): INDS = np.where(I == i)[0] A[INDS] = np.dot(A0[INDS], A1[i]) B[INDS] = np.dot(A0[INDS], B1[i]) + B0[INDS] return A, B @cached def _jacobian_inverse_transposed(self, codim): assert 0 <= codim < self.dim,\ f'Invalid Codimension (must be between 0 and {self.dim} but was {codim})' J = self.embeddings(codim)[0] if J.shape[-1] == J.shape[-2] == 2: JIT = inv_transposed_two_by_two(J) else: pinv = np.linalg.pinv JIT = np.array([pinv(j) for j in J]).swapaxes(1, 2) return JIT @cached def _integration_elements(self, codim): assert 0 <= codim <= self.dim,\ f'Invalid Codimension (must be between 0 and {self.dim} but was {codim})' if codim == self.dim: return np.ones(self.size(codim)) J = self.embeddings(codim)[0] JTJ = np.einsum('eji,ejk->eik', J, J) if JTJ.shape[1] == 1: D = JTJ.ravel() elif JTJ.shape[1] == 2: D = (JTJ[:, 0, 0] * JTJ[:, 1, 1] - JTJ[:, 1, 0] * JTJ[:, 0, 1]).ravel() else: def f(A): return np.linalg.det(A) D = np.array([f(j) for j in J]) return np.sqrt(D) @cached def _volumes(self, codim): assert 0 <= codim <= self.dim,\ f'Invalid Codimension (must be between 0 and {self.dim} but was {codim})' if codim == self.dim: return np.ones(self.size(self.dim)) return self.reference_element(codim).volume * self.integration_elements(codim) @cached def _volumes_inverse(self, codim): return np.reciprocal(self.volumes(codim)) @cached def _unit_outer_normals(self): JIT = self.jacobian_inverse_transposed(0) N = np.dot(JIT, self.reference_element(0).unit_outer_normals().T).swapaxes(1, 2) return N / np.apply_along_axis(np.linalg.norm, 2, N)[:, :, np.newaxis] @cached def _centers(self, codim): assert 0 <= codim <= self.dim,\ f'Invalid Codimension (must be between 0 and {self.dim} but was {codim})' A, B = self.embeddings(codim) C = self.reference_element(codim).center() return np.dot(A, C) + B @cached def _diameters(self, codim): assert 0 <= codim <= self.dim,\ f'Invalid Codimension (must be between 0 and {self.dim} but was {codim})' return np.reshape(self.reference_element(codim).mapped_diameter(self.embeddings(codim)[0]), (-1,)) @cached 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)) centers = self.centers(self.dim) for dim in range(self.dim): bbox[0, dim] = np.min(centers[:, dim]) bbox[1, dim] = np.max(centers[:, dim]) return bbox