# 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.analyticalproblems.domaindescriptions import KNOWN_BOUNDARY_TYPES
from pymor.core.base import abstractmethod
from pymor.core.cache import CacheableObject, cached
from pymor.discretizers.builtin.grids.defaultimpl import (ConformalTopologicalGridDefaultImplementations,
ReferenceElementDefaultImplementations,
AffineGridDefaultImplementations,)
[docs]class ReferenceElement(ReferenceElementDefaultImplementations, CacheableObject):
"""Defines a reference element.
All reference elements have the property that all subentities of a given codimension are of the
same type. I.e. a three-dimensional reference element cannot have triangles and rectangles as
faces at the same time.
Attributes
----------
dim
The dimension of the reference element
volume
The volume of the reference element
"""
dim = None
volume = None
cache_region = 'memory'
[docs] @abstractmethod
def size(self, codim):
"""Number of subentities of codimension `codim`."""
[docs] @abstractmethod
def subentities(self, codim, subentity_codim):
"""`subentities(c,sc)[i,j]` is, with respect to the indexing inside the
reference element, the index of the `j`-th codim-`subentity_codim`
subentity of the `i`-th codim-`codim` subentity of the reference element.
"""
pass
[docs] @abstractmethod
def subentity_embedding(self, subentity_codim):
"""Returns a tuple `(A, B)` which defines the embedding of the codim-`subentity_codim`
subentities into the reference element.
For `subentity_codim > 1', the embedding is by default given recursively via
`subentity_embedding(subentity_codim - 1)` and
`sub_reference_element(subentity_codim - 1).subentity_embedding(1)` choosing always
the superentity with smallest index.
"""
return self._subentity_embedding(subentity_codim)
[docs] @abstractmethod
def sub_reference_element(self, codim):
"""Returns the reference element of the codim-`codim` subentities."""
return self._sub_reference_element(codim)
[docs] def __call__(self, codim):
"""Returns the reference element of the codim-`codim` subentities."""
return self.sub_reference_element(codim)
[docs] @abstractmethod
def unit_outer_normals(self):
"""`retval[e]` is the unit outer-normal vector to the codim-1 subentity
with index `e`.
"""
pass
[docs] @abstractmethod
def center(self):
"""Coordinates of the barycenter."""
pass
[docs] @abstractmethod
def mapped_diameter(self, A):
"""The diameter of the reference element after transforming it with the
matrix `A` (vectorized).
"""
pass
[docs] @abstractmethod
def quadrature(self, order=None, npoints=None, quadrature_type='default'):
"""Returns tuple `(P, W)` where `P` is an array of quadrature points with
corresponding weights `W`.
The quadrature is of order `order` or has `npoints` integration points.
"""
pass
[docs] @abstractmethod
def quadrature_info(self):
"""Returns a tuple of dicts `(O, N)` where `O[quadrature_type]` is a list
of orders which are implemented for `quadrature_type` and `N[quadrature_type]`
is a list of the corresponding numbers of integration points.
"""
pass
def quadrature_types(self):
o, _ = self.quadrature_info()
return frozenset(o.keys())
[docs]class AffineGrid(AffineGridDefaultImplementations, ConformalTopologicalGrid):
"""Topological grid with geometry where each codim-0 entity is affinely mapped to the same |ReferenceElement|.
The grid is completely determined via the subentity relation given by
:meth:`~AffineGrid.subentities` and the embeddings given by
:meth:`~AffineGrid.embeddings`. In addition,
only :meth:`~ConformalTopologicalGrid.size` and :meth:`~AffineGrid.reference_element`
have to be implemented. Cached default implementations for all other methods are
provided by :class:`~pymor.discretizers.builtin.grids.defaultimpl.AffineGridDefaultImplementations`.
"""
[docs] @abstractmethod
def reference_element(self, codim):
"""The |ReferenceElement| of the codim-`codim` entities."""
pass
[docs] @abstractmethod
def subentities(self, codim, subentity_codim):
"""`retval[e,s]` is the global index of the `s`-th codim-`subentity_codim`
subentity of the codim-`codim` entity with global index `e`.
The ordering of `subentities(0, subentity_codim)[e]` has to correspond, w.r.t.
the embedding of `e`, to the local ordering inside the reference element.
For `codim > 0`, we provide a default implementation by calculating the
subentities of `e` as follows:
1. Find the `codim-1` parent entity `e_0` of `e` with minimal global index
2. Lookup the local indices of the subentities of `e` inside `e_0` using the reference element.
3. Map these local indices to global indices using `subentities(codim - 1, subentity_codim)`.
This procedures assures that `subentities(codim, subentity_codim)[e]`
has the right ordering w.r.t. the embedding determined by `e_0`, which
agrees with what is returned by `embeddings(codim)`
"""
return self._subentities(codim, subentity_codim)
[docs] @abstractmethod
def embeddings(self, codim):
"""Returns tuple `(A, B)` where `A[e]` and `B[e]` are the linear part
and the translation part of the map from the reference element of `e`
to `e`.
For `codim > 0`, we provide a default implementation by
taking the embedding of the codim-1 parent entity `e_0` of `e` with
lowest global index and composing it with the subentity_embedding of `e`
into `e_0` determined by the reference element.
"""
return self._embeddings(codim)
[docs] def jacobian_inverse_transposed(self, codim):
"""`retval[e]` is the transposed (pseudo-)inverse of the Jacobian of `embeddings(codim)[e]`.
"""
return self._jacobian_inverse_transposed(codim)
[docs] def integration_elements(self, codim):
"""`retval[e]` is given as `sqrt(det(A^T*A))`, where `A = embeddings(codim)[0][e]`."""
return self._integration_elements(codim)
[docs] def volumes(self, codim):
"""`retval[e]` is the (dim-`codim`)-dimensional volume of the codim-`codim` entity with global index `e`."""
return self._volumes(codim)
[docs] def volumes_inverse(self, codim):
"""`retval[e] = 1 / volumes(codim)[e]`."""
return self._volumes_inverse(codim)
[docs] def unit_outer_normals(self):
"""`retval[e,i]` is the unit outer normal to the i-th codim-1 subentity
of the codim-0 entitiy with global index `e`.
"""
return self._unit_outer_normals()
[docs] def centers(self, codim):
"""`retval[e]` is the barycenter of the codim-`codim` entity with global index `e`."""
return self._centers(codim)
[docs] def diameters(self, codim):
"""`retval[e]` is the diameter of the codim-`codim` entity with global index `e`."""
return self._diameters(codim)
[docs] def quadrature_points(self, codim, order=None, npoints=None, quadrature_type='default'):
"""`retval[e]` is an array of quadrature points in global coordinates
for the codim-`codim` entity with global index `e`.
The quadrature is of order `order` or has `npoints` integration points. To
integrate a function `f` over `e` one has to form ::
np.dot(f(quadrature_points(codim, order)[e]), reference_element(codim).quadrature(order)[1]) *
integration_elements(codim)[e]. # NOQA
"""
return self._quadrature_points(codim, order, npoints, quadrature_type)
[docs] def bounding_box(self):
"""returns a `(2, dim)`-shaped array containing lower/upper bounding box coordinates."""
return self._bounding_box()
[docs]class AffineGridWithOrthogonalCenters(AffineGrid):
"""|AffineGrid| with an additional `orthogonal_centers` method."""
[docs] @abstractmethod
def orthogonal_centers(self):
"""`retval[e]` is a point inside the codim-0 entity with global index `e`
such that the line segment from `retval[e]` to `retval[e2]` is always
orthogonal to the codim-1 entity shared by the codim-0 entities with global
index `e` and `e2`.
(This is mainly useful for gradient approximation in finite volume schemes.)
"""
pass
[docs]class BoundaryInfo(CacheableObject):
"""Provides boundary types for the boundaries of a given |ConformalTopologicalGrid|.
For every boundary type and codimension a mask is provided, marking grid entities
of the respective type and codimension by their global index.
Attributes
----------
boundary_types
set of all boundary types the grid has.
"""
boundary_types = frozenset()
cache_region = 'memory'
[docs] def mask(self, boundary_type, codim):
"""retval[i] is `True` if the codim-`codim` entity of global index `i` is
associated to the boundary type `boundary_type`.
"""
raise ValueError(f'Has no boundary_type "{boundary_type}"')
[docs] def unique_boundary_type_mask(self, codim):
"""retval[i] is `True` if the codim-`codim` entity of global index `i` is
associated to one and only one boundary type.
"""
return np.less_equal(sum(self.mask(bt, codim=codim).astype(np.int) for bt in self.boundary_types), 1)
[docs] def no_boundary_type_mask(self, codim):
"""retval[i] is `True` if the codim-`codim` entity of global index `i` is
associated to no boundary type.
"""
return np.equal(sum(self.mask(bt, codim=codim).astype(np.int) for bt in self.boundary_types), 0)
def check_boundary_types(self, assert_unique_type=(1,), assert_some_type=()):
for bt in self.boundary_types:
if bt not in KNOWN_BOUNDARY_TYPES:
self.logger.warning(f'Unknown boundary type: {bt}')
if assert_unique_type:
for codim in assert_unique_type:
assert np.all(self.unique_boundary_type_mask(codim))
if assert_some_type:
for codim in assert_some_type:
assert not np.any(self.no_boundary_type_mask(codim))
@property
def has_dirichlet(self):
return 'dirichlet' in self.boundary_types
@property
def has_neumann(self):
return 'neumann' in self.boundary_types
@property
def has_robin(self):
return 'robin' in self.boundary_types
def dirichlet_mask(self, codim):
return self.mask('dirichlet', codim)
def neumann_mask(self, codim):
return self.mask('neumann', codim)
def robin_mask(self, codim):
return self.mask('robin', codim)
@cached
def _boundaries(self, boundary_type, codim):
return np.where(self.mask(boundary_type, codim))[0].astype('int32')
def boundaries(self, boundary_type, codim):
return self._boundaries(boundary_type, codim)
def dirichlet_boundaries(self, codim):
return self._boundaries('dirichlet', codim)
def neumann_boundaries(self, codim):
return self._boundaries('neumann', codim)
def robin_boundaries(self, codim):
return self._boundaries('robin', codim)