Source code for sfepy.discrete.dg.poly_spaces

# -*- coding: utf-8 -*-
from abc import abstractmethod
from functools import reduce
from operator import mul

import numpy as nm

from sfepy.base.ioutils import InDir
from sfepy.base.base import Struct
from sfepy.discrete import PolySpace
from scipy.special import jacobi as scp_jacobi
from scipy.special import eval_jacobi as scp_eval_jacobi


[docs] def iter_by_order(order, dim, extended=False): """Iterates over all combinations of basis functions indexes needed to create multidimensional basis in a way that creates hierarchical basis Parameters ---------- order : int desired order of multidimensional basis dim : int dimension of the basis extended : bool iterate over extended tensor product basis (Default value = False) Yields ------ idx : tuple containing basis function indexes, used in ``_combine_polyvals`` and ``_combine_polyvals_der`` """ # nth(iter(map(lambda x: x + (order - reduce(add,x),)), range(order)), dim) # nth(dim, iterate(map(lambda x: x + (order - reduce(add,x),)), # map(tuple, range(order)))) # nth(2, iterate(map(lambda x: x + (order - reduce(add,x),)), # map(lambda x: (x,), range(order)))) porder = order + 1 if dim == 1: for i in range(porder): yield i, return elif dim == 2: for k in range(porder): for i in range(k + 1): yield i, k - i if not extended: return for k in range(1, porder): for i in range(1, porder): if i + k <= porder - 1: continue yield i, k elif dim == 3: for k in range(porder): for j in range(k + 1): for i in range(j + 1): yield i, j - i, k - j return
[docs] def get_n_el_nod(order, dim, extended=False): r"""Number of nodes per element for discontinuous legendre basis, i.e. number of iterations yielded by iter_by_order When extended is False .. math:: N_p = \frac{(n + 1) \cdot (n + 2) \cdot ... \cdot (n + d)}{d!} where `n` is the order and `d` the dimension. When extended is True .. math:: N_p = (n + 1) ^ d where `n` is the order and `d` the dimension. Parameters ---------- order : int desired order of multidimensional basis dim : int dimension of the basis extended : bool iterate over extended tensor product basis (Default value = False) Returns ------- n_el_nod : int number of basis functions in basis """ if extended: return (order + 1) ** dim return int(reduce(mul, map(lambda i: order + i + 1, range(dim))) / reduce(mul, range(1, dim + 1)))
[docs] class LegendrePolySpace(PolySpace): """Legendre hierarchical polynomials basis, over [0, 1] domain.""" def __init__(self, name, geometry, order, extended): """ Parameters ---------- name geometry order : int approximation order, 0 for constant functions, 1 for linear etc. extended : bool for extended tensor product space """ PolySpace.__init__(self, name, geometry, order) self.extended = extended # only tensor product polyspace is extended self.n_v = geometry.n_vertex, self.dim = geometry.dim self.n_nod = get_n_el_nod(self.order, self.dim, self.extended) self.coefM = None self.expoM = None def _eval_basis(self, coors, diff=0, ori=None, suppress_errors=False, eps=1e-15): """Calls _combine_polyvals or _combine_polyvals_diff to build multidimensional basis implemented in subclasses to get basis for different geometries expects coors to be in shape (..., dim), Returns values of the basis in shape (coors.shape[:-1], 1, n_el_nod) or values of the gradient in shape (1, coors.shape[:-1], dim, n_el_nod) Parameters ---------- coors : array_like ori : unused we do not need ori, because the basis is discontinous across elements this is treated in DGField(Default value = None) suppress_errors : has no effect diff : int or truthy (Default value = 0) eps : unused (Default value = 1e-15) Returns ------- values : ndarray """ coors = 2 * coors - 1 # transofrm from [0, 1] to [-1, 1] porder = self.order + 1 if diff: diff = int(diff) values = nm.zeros((1,) + coors.shape[:-1] + # (1,) for dummy axis used throughout sfepy (self.dim,) * diff + # (dim,) * diff order is shape of derivation tensor, # we support only first derivative (self.n_nod,)) polyvals = nm.zeros(coors.shape + (porder,) + (diff + 1,)) # diff + 1 is number of values of one dimensional base polyvals[..., 0] = self.legendreP(coors) polyvals[..., 1] = self.gradlegendreP(coors) for m, idx in enumerate( iter_by_order(self.order, self.dim, self.extended)): for d in range(self.dim): values[..., d, m] = 2 * self._combine_polyvals_diff(coors, polyvals, d, idx) # 2 is due to transformation from [0,1] to [-1,1] else: values = nm.zeros(coors.shape[:-1] + (1, self.n_nod,)) # 1, because no matter the dimension functions have only one value polyvals = self.legendreP(coors) for m, idx in enumerate( iter_by_order(self.order, self.dim, self.extended)): values[..., 0, m] = self._combine_polyvals(coors, polyvals, idx) return values
[docs] def get_interpol_scheme(self): r"""For dim > 1 returns F and P matrices according to gmsh basis specification [1]_: Let us assume that the approximation of the view's value over an element is written as a linear combination of d basis functions :math:`f_i, i=0, ..., n-1` (the coefficients being stored in list-of-values). Defining .. math :: f_i = \sum\limits_{j=0}^{d-1} F_{ij}\cdot p_j, with :math:`p_j(u, v, w) = u^{P_j^{(0)}}\cdot v^{P_j^{(1)}}\cdot w^{P_j^{(2)}}` (`u`, `v` and `w` being the coordinates in the element's parameter space), then val-coef-matrix denotes the n x n matrix F and val-exp-matrix denotes the n x 3 matrix P where n is number of basis functions as calculated by ``get_n_el_nod``. Expects matrices to be saved in atributes coefM and expoM! .. [1] Remacle, J.-F., Chevaugeon, N., Marchandise, E., & Geuzaine, C. (2007). Efficient visualization of high-order finite elements. International Journal for Numerical Methods in Engineering, 69(4), 750-771. https://doi.org/10.1002/nme.1787 Returns ------- interp_scheme_struct : Struct Struct with name of the scheme, geometry desc and P and F """ interp_scheme_struct = None if self.dim > 1: interp_scheme_struct = Struct(name=self.name + "_interpol_scheme", F=self.coefM, P=self.expoM, desc=self.geometry.name) return interp_scheme_struct
@abstractmethod def _combine_polyvals(self, coors, polyvals, idx): """ Combines Legendre or Jacobi polynomials to get muldidimensional basis values according to element topolgy Parameters ---------- coors : array_like coordinates of evaluation points polyvals : array_like values of legendre polynomials precomputed in _eval_basis idx : tuple function index, for tensor-product correspond to orders of polynomials in variables Returns ------- values : ndarray """ @abstractmethod def _combine_polyvals_diff(self, coors, polyvals, der, idx): """Combines Legendre or Jacobi polynomials to get muldidimensional basis derivative values according to element topolgy. Parameters ---------- coors : array_like coordinates of evaluation points polyvals : array_like values of legendre polynomials precomputed in _eval_basis der : int derivative variable idx : tuple function index, for tensor-product correspond to orders of polynomials in variables Returns ------- values : ndarray """ # --------------------- # # 1D legendre polyspace # # --------------------- #
[docs] def legendreP(self, coors): """ Parameters ---------- coors : array_like coordinates, preferably in interval [-1, 1] for which this basis is intented Returns ------- values : ndarray values at coors of all the legendre polynomials up to self.order """ return self.jacobiP(coors, alpha=0, beta=0)
[docs] def gradlegendreP(self, coors, diff=1): """ Parameters ---------- diff : int default 1 coors : array_like coordinates, preferably in interval [-1, 1] for which this basis is intented Returns ------- values : ndarray values at coors of all the legendre polynomials up to self.order """ return self.gradjacobiP(coors, 0, 0, diff=diff)
""" Explict legendre polynomials up to order 5 """ legendre_funs = [lambda x: 0 * x + 1, # we need constant preserving shape and type of x lambda x: 2 * x - 1, lambda x: (6 * x ** 2 - 6 * x + 1), lambda x: (20 * x ** 3 - 30 * x ** 2 + 12 * x - 1), lambda x: (70 * x ** 4 - 140 * x ** 3 + 90 * x ** 2 - 20 * x + 1), lambda x: (252 * x ** 5 - 630 * x ** 4 + 560 * x ** 3 - 210 * x ** 2 + 30 * x - 1) ]
[docs] def get_nth_fun(self, n): """Uses shifted Legendre polynomials formula on interval [0, 1]. Convenience function for testing Parameters ---------- n : int Returns ------- fun : callable n-th function of the legendre basis """ if n < 6: return self.legendre_funs[n] else: from scipy.special import comb as comb def fun(x): val = 0.0 for k in range(n + 1): val = val + comb(n, k) * comb(n + k, k) * (-x) ** k val *= -1 if n % 2 else 1 return val return fun
[docs] def get_nth_fun_der(self, n, diff=1): """Returns diff derivative of nth function. Uses shifted legendre polynomials formula on interval [0, 1]. Useful for testing. Parameters ---------- n : int diff : int (Default value = 1) Returns ------- fun : callable derivative of n-th function of the 1D legendre basis """ def dfun(x): """ Parameters ---------- x : Returns ------- """ from scipy.special import comb as comb, factorial val = x[:] * 0.0 for k in range(diff, n + 1): val += comb(n, k) * comb(n + k, k) * factorial(k) / \ factorial(k - diff) * (-x) ** (k - diff) val *= -1 if (n - diff) % 2 else 1 return val return dfun
# --------------------- # # --------------------- # # 1D jacobi polyspace # # --------------------- #
[docs] def jacobiP(self, coors, alpha, beta): """Values of the jacobi polynomials on interval [-1, 1] up to self.order + 1 at coors Parameters ---------- coors : array_like beta : float alpha : float Returns ------- values : ndarray output shape is shape(coor) + (self.order + 1,) """ if not isinstance(coors, nm.ndarray): sh = (1,) else: sh = nm.shape(coors) values = nm.ones(sh + (self.order + 1,)) for i in range(self.order + 1): values[..., i] = scp_eval_jacobi(i, alpha, beta, coors) # for some reason eval_jacobi consumes last dimension if it is one, # when called with order array return values
[docs] def gradjacobiP(self, coors, alpha, beta, diff=1): """diff derivative of the jacobi polynomials on interval [-1, 1] up to self.order + 1 at coors Parameters ---------- coors : alpha : float beta : float diff : int (Default value = 1) Returns ------- values : ndarray output shape is shape(coor) + (self.order + 1,) """ if isinstance(coors, (int, float)): sh = (1,) else: sh = nm.shape(coors) values = nm.zeros(sh + (self.order + 1,)) for i in range(self.order + 1): jacob_poly = scp_jacobi(i, alpha, beta) values[..., i] = jacob_poly.deriv(m=diff)(coors) # Warning # Computing values of high-order polynomials (around order > 20) # using polynomial coefficients is numerically unstable. To evaluate # polynomial values, the eval_* functions should be used instead. # On that note jacob_poly.deriv seems to use stable version. return values
# --------------------- #
[docs] class LegendreTensorProductPolySpace(LegendrePolySpace): name = "legendre_tensor_product" def __init__(self, name, geometry, order, ): super().__init__(name, geometry, order, extended=True) self.n_nod = get_n_el_nod(self.order, self.dim, self.extended) if self.dim > 1: self.coefM, self.expoM = self.build_interpol_scheme() def _combine_polyvals(self, coors, polyvals, idx): return nm.prod(polyvals[..., range(len(idx)), idx], axis=-1) def _combine_polyvals_diff(self, coors, polyvals, dvar, idx): dimz = range(len(idx)) derz = nm.zeros(len(idx), dtype=nm.int32) derz[dvar] = 1 return nm.prod(polyvals[..., dimz, idx, derz], axis=-1)
[docs] def build_interpol_scheme(self): """Builds F and P matrices returned by self.get_interpol_scheme. Note that this function returns coeficients according to gmsh parametrization of Quadrangle i.e. [-1, 1] x [-1, 1] and hence the form of basis function is not the same as exhibited by the LegendreTensorProductPolySpace object which acts on parametrization [0, 1] x [0, 1]. Returns ------- F : ndarray coefficient matrix P : ndarray exponent matrix """ P = nm.zeros((self.n_nod, 3), dtype=nm.int32) for m, idx in enumerate( iter_by_order(self.order, self.dim, self.extended)): P[m, :self.dim] = idx F = nm.zeros((self.n_nod, self.n_nod)) for m, idx in enumerate( iter_by_order(self.order, self.dim, self.extended)): xcoefs = list(scp_jacobi(idx[0], 0, 0).coef)[::-1] xcoefs = nm.array(xcoefs + [0] * (self.order + 1 - len(xcoefs))) ycoefs = list(scp_jacobi(idx[1], 0, 0).coef)[::-1] ycoefs = nm.array(ycoefs + [0] * (self.order + 1 - len(ycoefs))) coef_mat = nm.outer(xcoefs, ycoefs) F[m, :] = [coef_mat[idx] for idx in iter_by_order(self.order, self.dim, self.extended)] return F, P
[docs] class LegendreSimplexPolySpace(LegendrePolySpace): name = "legendre_simplex" def __init__(self, name, geometry, order, extended=False): super().__init__(name, geometry, order, extended) if self.dim == 1: return if order <= 5: self.coefM = simplex_coefM5[:self.n_nod, :self.n_nod] self.expoM = simplex_expoM5[:self.n_nod, :] return indir = InDir(__file__) try: self.coefM = nm.loadtxt( indir("legendre2D_simplex_coefs.txt") )[:self.n_nod, :self.n_nod] self.expoM = nm.loadtxt( indir("legendre2D_simplex_expos.txt") )[:self.n_nod, :] except IOError as e: raise IOError( ("File {} not found, run gen_legendre_simplex_base.py" + " to generate it. This is needed approx order > 5.") .format(e.args[0])) def _combine_polyvals(self, coors, polyvals, idx): from scipy.special import eval_jacobi if len(idx) == 1: # 1D return nm.prod(polyvals[..., range(len(idx)), idx], axis=-1) elif len(idx) == 2: # 2D r = coors[..., 0] s = coors[..., 1] a = 2 * (1 + r) / (1 - s) - 1 a[nm.isnan(a)] = -1. b = s return eval_jacobi(idx[0], 0, 0, a) \ * eval_jacobi(idx[1],2 * idx[0] + 1, 0,b) * (1 - b) ** idx[0] elif len(idx) == 3: # 3D r = coors[..., 0] s = coors[..., 1] t = coors[..., 2] a = -2 * (1 + r) / (s + t) - 1 b = 2 * (1 + s) / (1 - t) - t c = t a[nm.isnan(a)] = -1. b[nm.isnan(b)] = -1. return eval_jacobi(idx[0], 0, 0, a) * \ eval_jacobi(idx[1], 2 * idx[0] + 1, 0, 0, b) * \ eval_jacobi(idx[2], 2 * idx[0] + 2 * idx[1] + 2, 0, c) * \ (1 - c) ** (idx[0] + idx[1]) def _combine_polyvals_diff(self, coors, polyvals, dvar, idx): if len(idx) == 1: # 1D dimz = range(len(idx)) derz = nm.zeros(len(idx), dtype=nm.int32) derz[dvar] = 1 return nm.prod(polyvals[..., dimz, idx, derz], axis=-1) elif len(idx) == 2: # 2D r = coors[..., 0] s = coors[..., 1] a = 2 * (1 + r) / (1 - s) - 1 b = s a[nm.isnan(a)] = -1. di = idx[0] dj = idx[1] fa = self.jacobiP(a, 0, 0)[..., di] dfa = self.gradjacobiP(a, 0, 0)[..., di] gb = self.jacobiP(b, 2 * di + 1, 0)[..., dj] dgb = self.gradjacobiP(b, 2 * di + 1, 0)[..., dj] if dvar == 0: # d/dr # r - derivative # d / dr = da / dr * d / da + db / dr * d / db # = (2 / (1 - s)) d / da = (2 / (1 - b)) d / da dmodedr = dfa * gb dmodedr = dmodedr * ( (0.5 * (1 - b)) ** (di - 1)) if di > 0 else dmodedr return 2 ** di * dmodedr elif dvar == 1: # d/ds # s - derivative # d / ds = ((1 + a) / 2) / ((1 - b) / 2) d / da + d / db dmodeds = dfa * (gb * (0.5 * (1 + a))) dmodeds = dmodeds * ( (0.5 * (1 - b)) ** (di - 1)) if di > 0 else dmodeds tmp = dgb * ((0.5 * (1 - b)) ** di) tmp = tmp - 0.5 * di * gb * ( (0.5 * (1 - b)) ** (di - 1)) if di > 0 else tmp dmodeds = dmodeds + fa * tmp return 2 ** di * dmodeds elif len(idx) == 3: # 3D # UNTESTED! r = coors[..., 0] s = coors[..., 1] t = coors[..., 2] a = -2 * (1 + r) / (s + t) - 1 b = 2 * (1 + s) / (1 - t) - t c = t a[nm.isnan(a)] = -1. b[nm.isnan(b)] = -1. di = idx[0] dj = idx[1] dk = idx[2] fa = self.jacobiP(a, 0, 0)[..., di] dfa = self.gradjacobiP(a, 0, 0)[..., di] gb = self.jacobiP(b, 2 * di + 1, 0)[..., dj] dgb = self.gradjacobiP(b, 2 * di + 1, 0)[..., dj] hc = self.jacobiP(c, 2 * (di + dj) + 2, 0, dk) dhc = self.gradjacobiP(c, 2 * (di + dj) + 2, 0, dk) V3Dr = dfa * (gb * hc) V3Dr = V3Dr * ((0.5 * (1 - b)) ** (di - 1)) if di > 0 else V3Dr V3Dr = V3Dr * ((0.5 * (1 - c)) ** ( di + dj - 1)) if di + dj > 0 else V3Dr tmp = dgb * ((0.5 * (1 - b)) ** di) tmp = tmp + (-0.5 * di) * ( gb * (0.5 * (1 - b)) ** (di - 1)) if di > 0 else tmp tmp = tmp * ((0.5 * (1 - c)) ** ( di + dj - 1)) if di + dj > 0 else tmp tmp = fa * (tmp * hc) if dvar == 0: return V3Dr * (2 ** (2 * di + dj + 1)) elif dvar == 1: V3Ds = 0.5 * (1 + a) * V3Dr V3Ds = V3Ds + tmp return V3Ds * (2 ** (2 * di + dj + 1)) elif dvar == 2: V3Dt = 0.5 * (1 + a) * V3Dr + 0.5 * (1 + b) * tmp tmp = dhc * ((0.5 * (1 - c)) ** (di + dj)) tmp = tmp - 0.5 * (di + dj) * (hc * ((0.5 * (1 - c)) ** ( di + dj - 1))) if di + dj > 0 else tmp tmp = fa * (gb * tmp) tmp = tmp * ((0.5 * (1 - b)) ** di) V3Dt = V3Dt + tmp return V3Dt * (2 ** (2 * di + dj + 1))
simplex_coefM5 = nm.array([ 1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, -1,3,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, -2,2,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 1,-8,0,10,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 2,-12,-4,10,20,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 4,-8,-24,4,24,24,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0, -1,15,0,-45,0,0,35,0,0,0,0,0,0,0,0,0,0,0,0,0,0, -2,26,4,-66,-48,0,42,84,0,0,0,0,0,0,0,0,0,0,0,0,0, -4,36,24,-60,-192,-24,28,168,168,0,0,0,0,0,0,0,0,0,0,0,0, -8,24,96,-24,-192,-240,8,96,240,160,0,0,0,0,0,0,0,0,0,0,0, 1,-24,0,126,0,0,-224,0,0,0,126,0,0,0,0,0,0,0,0,0,0, 2,-44,-4,210,84,0,-336,-336,0,0,168,336,0,0,0,0,0,0,0,0,0, 4,-72,-24,276,408,24,-352,-1248,-384,0,144,864,864,0,0,0,0,0,0,0,0, 8,-96,-96,240,1056,240,-224,-1824,-2400,-160,72,864,2160,1440,0,0,0,0,0,0,0, 16,-64,-320,96,960,1440,-64,-960,-2880,-2240,16,320,1440,2240,1120,0,0,0,0,0,0, -1,35,0,-280,0,0,840,0,0,0,-1050,0,0,0,0,462,0,0,0,0,0, -2,66,4,-496,-128,0,1392,864,0,0,-1620,-1920,0,0,0,660,1320,0,0,0,0, -4,116,24,-760,-672,-24,1848,3888,648,0,-1860,-7200,-3240,0,0,660,3960,3960,0,0,0, -8,184,96,-944,-2112,-240,1808,9216,5040,160,-1480,-12480,-18000,-3200,0,440,5280,13200,8800,0,0, -16,240,320,-800,-4480,-1440,1120,11520,18720,2240,-720,-10880,-33120,-26880,-1120,176,3520,15840,24640,12320,0, -32,160,960,-320,-3840,-6720,320,5760,20160,17920,-160,-3840,-20160,-35840,-20160,32,960,6720,17920,20160,8064], dtype=nm.int64).reshape(21, 21) simplex_expoM5 = nm.array([ 0,0,0, 0,1,0, 1,0,0, 0,2,0, 1,1,0, 2,0,0, 0,3,0, 1,2,0, 2,1,0, 3,0,0, 0,4,0, 1,3,0, 2,2,0, 3,1,0, 4,0,0, 0,5,0, 1,4,0, 2,3,0, 3,2,0, 4,1,0, 5,0,0], dtype=nm.int64).reshape(21, 3)