from __future__ import absolute_import
import numpy as nm
from sfepy.base.base import load_classes, Struct
from sfepy import get_paths
def _get_table():
if PolySpace._all is None:
ps_files = get_paths('sfepy/discrete/fem/poly_spaces.py')
ps_files += get_paths('sfepy/discrete/dg/poly_spaces.py')
PolySpace._all = load_classes(ps_files, [PolySpace],
ignore_errors=True,
name_attr='name')
return PolySpace._all
[docs]
def register_poly_space(cls):
table = _get_table()
table[cls.name] = cls
[docs]
class PolySpace(Struct):
"""Abstract polynomial space class."""
_all = None
keys = {
(0, 1) : 'simplex',
(1, 2) : 'simplex',
(2, 3) : 'simplex',
(3, 4) : 'simplex',
(3, 6) : 'wedge',
(2, 4) : 'tensor_product',
(3, 8) : 'tensor_product',
}
[docs]
@staticmethod
def any_from_args(name, geometry, order, basis='lagrange',
force_bubble=False):
"""
Construct a particular polynomial space classes according to the
arguments passed in.
"""
if name is None:
name = PolySpace.suggest_name(geometry, order, basis, force_bubble)
table = _get_table()
key = '%s_%s' % (basis, PolySpace.keys[(geometry.dim,
geometry.n_vertex)])
if (geometry.name == '1_2') and (key not in table):
key = '%s_%s' % (basis, 'tensor_product')
if force_bubble:
key += '_bubble'
return table[key](name, geometry, order)
[docs]
@staticmethod
def suggest_name(geometry, order, basis='lagrange',
force_bubble=False):
"""
Suggest the polynomial space name given its constructor parameters.
"""
aux = geometry.get_interpolation_name()[:-1]
if force_bubble:
return aux + ('%dB' % order)
else:
return aux + ('%d' % order)
def __init__(self, name, geometry, order):
self.name = name
self.geometry = geometry
self.order = order
self.bbox = nm.vstack((geometry.coors.min(0), geometry.coors.max(0)))
[docs]
def eval_basis(self, coors, diff=0, ori=None, force_axis=False,
transform=None, suppress_errors=False, eps=1e-15):
"""
Evaluate the basis or its first or second derivatives in points given
by coordinates. The real work is done in _eval_basis() implemented in
subclasses.
Note that the second derivative code is a work-in-progress and only
`coors` and `transform` arguments are used.
Parameters
----------
coors : array_like
The coordinates of points where the basis is evaluated. See Notes.
diff : 0, 1 or 2
If nonzero, return the given derivative.
ori : array_like, optional
Optional orientation of element facets for per element basis.
force_axis : bool
If True, force the resulting array shape to have one more axis even
when `ori` is None.
transform : array_like, optional
The basis transform array.
suppress_errors : bool
If True, do not report points outside the reference domain.
eps : float
Accuracy for comparing coordinates.
Returns
-------
basis : array
The basis (shape (n_coor, 1, n_basis)) or its first derivative
(shape (n_coor, dim, n_basis)) or its second derivative (shape
(n_coor, dim, dim, n_basis)) evaluated in the given points. An
additional axis is pre-pended of length n_cell, if `ori` is given,
or of length 1, if `force_axis` is True.
Notes
-----
If coors.ndim == 3, several point sets are assumed, with equal number
of points in each of them. This is the case, for example, of the
values of the volume basis functions on the element facets. The indexing
(of bf_b(g)) is then (ifa,iqp,:,n_ep), so that the facet can be set in
C using FMF_SetCell.
"""
coors = nm.asarray(coors)
if not coors.ndim in (2, 3):
raise ValueError('coordinates must have 2 or 3 dimensions! (%d)'
% coors.ndim)
if (coors.shape[-1] != self.geometry.dim) and (self.geometry.dim > 0):
raise ValueError('PolySpace geometry dimension %d does not agree'
' with quadrature coordinates dimension %d!'
% (self.geometry.dim, coors.shape[-1]))
if (coors.ndim == 2):
basis = self._eval_basis(coors, diff=diff, ori=ori,
suppress_errors=suppress_errors,
eps=eps)
if (basis.ndim == 3) and force_axis:
basis = basis[None, ...]
if not basis.flags['C_CONTIGUOUS']:
basis = nm.ascontiguousarray(basis)
else: # Several point sets.
if diff:
bdim = self.geometry.dim
else:
bdim = 1
basis = nm.empty((coors.shape[0], coors.shape[1],
bdim, self.n_nod), dtype=nm.float64)
for ii, _coors in enumerate(coors):
basis[ii] = self._eval_basis(_coors, diff=diff, ori=ori,
suppress_errors=suppress_errors,
eps=eps)
if transform is not None:
basis = transform_basis(transform, basis)
return basis