diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index 4f0a3a01482af78dacbfe0635ce45249a610f930..e1da810e632822e67e167470e1f02e87538b5b67 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -162,11 +162,12 @@ class ExpansionBase(object):
class ExpansionTermsWrangler(object):
- init_arg_names = ("order", "dim")
+ init_arg_names = ("order", "dim", "max_mi")
- def __init__(self, order, dim):
+ def __init__(self, order, dim, max_mi=None):
self.order = order
self.dim = dim
+ self.max_mi = max_mi
def get_coefficient_identifiers(self):
raise NotImplementedError
@@ -189,6 +190,16 @@ class ExpansionTermsWrangler(object):
as gnitstam)
res = sorted(gnitstam(self.order, self.dim), key=sum)
+
+ def filter_tuple(tup):
+ if self.max_mi is None:
+ return True
+ for a, b in zip(tup, self.max_mi):
+ if a > b:
+ return False
+ return True
+
+ res = list(filter(filter_tuple, res))
return res
def copy(self, **kwargs):
@@ -196,8 +207,8 @@ class ExpansionTermsWrangler(object):
(name, getattr(self, name))
for name in self.init_arg_names)
- new_kwargs["order"] = kwargs.pop("order", self.order)
- new_kwargs["dim"] = kwargs.pop("dim", self.dim)
+ for name in self.init_arg_names:
+ new_kwargs[name] = kwargs.pop(name, getattr(self, name))
if kwargs:
raise TypeError("unexpected keyword arguments '%s'"
@@ -261,9 +272,9 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler):
.. automethod:: get_reduced_coeffs
"""
- init_arg_names = ("order", "dim", "deriv_multiplier")
+ init_arg_names = ("order", "dim", "deriv_multiplier", "max_mi")
- def __init__(self, order, dim, deriv_multiplier):
+ def __init__(self, order, dim, deriv_multiplier, max_mi=None):
r"""
:param order: order of the expansion
:param dim: number of dimensions
@@ -274,7 +285,7 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler):
that representation of the PDE is :math:`\text{coeff} /
{\text{deriv\_multiplier}^{|\nu|}}`.
"""
- super(LinearRecurrenceBasedExpansionTermsWrangler, self).__init__(order, dim)
+ super(LinearRecurrenceBasedExpansionTermsWrangler, self).__init__(order, dim, max_mi)
self.deriv_multiplier = deriv_multiplier
def get_coefficient_identifiers(self):
@@ -354,7 +365,6 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler):
pde_dict = self.get_pde_dict()
pde_mat = []
-
mis = self.get_full_coefficient_identifiers()
coeff_ident_enumerate_dict = dict((tuple(mi), i) for
(i, mi) in enumerate(mis))
@@ -385,11 +395,10 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler):
Then, :math:`S^T = N * N_{[r, :]}^{-1}` which implies,
:math:`S = N_{[r, :]}^{-T} N^T = N_{[r, :]}^{-T} N^T`.
"""
- pde_mat = np.array(pde_mat, dtype=np.float64)
- n = nullspace(pde_mat, atol=tol)
+ n = nullspace(pde_mat)
idx = self.get_reduced_coeffs()
assert len(idx) >= n.shape[1]
- s = np.linalg.solve(n.T[:, idx], n.T)
+ s = n.T[:, idx].solve(n.T)
stored_identifiers = [mis[i] for i in idx]
else:
s = np.eye(len(mis))
@@ -400,8 +409,8 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler):
coeff_matrix = defaultdict(list)
for i in range(s.shape[0]):
for j in range(s.shape[1]):
- if np.abs(s[i][j]) > tol:
- coeff_matrix[j].append((i, s[i][j]))
+ if np.abs(s[i, j]) > tol:
+ coeff_matrix[j].append((i, s[i, j]))
plog.done()
@@ -453,10 +462,10 @@ class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler):
class LaplaceExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler):
- init_arg_names = ("order", "dim")
+ init_arg_names = ("order", "dim", "max_mi")
- def __init__(self, order, dim):
- super(LaplaceExpansionTermsWrangler, self).__init__(order, dim, 1)
+ def __init__(self, order, dim, max_mi=None):
+ super(LaplaceExpansionTermsWrangler, self).__init__(order, dim, 1, max_mi)
def get_pde_dict(self):
pde_dict = {}
@@ -478,13 +487,13 @@ class LaplaceExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler)
class HelmholtzExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler):
- init_arg_names = ("order", "dim", "helmholtz_k_name")
+ init_arg_names = ("order", "dim", "helmholtz_k_name", "max_mi")
- def __init__(self, order, dim, helmholtz_k_name):
+ def __init__(self, order, dim, helmholtz_k_name, max_mi=None):
self.helmholtz_k_name = helmholtz_k_name
multiplier = sym.Symbol(helmholtz_k_name)
- super(HelmholtzExpansionTermsWrangler, self).__init__(order, dim, multiplier)
+ super(HelmholtzExpansionTermsWrangler, self).__init__(order, dim, multiplier, max_mi)
def get_pde_dict(self, **kwargs):
pde_dict = {}
diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index 11b1368a43bd8c57ad51065f84bea3b7dc2a5ca5..cf94ebd674201178047cfee86acca0d74512c1ec 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -167,11 +167,18 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
from sumpy.tools import add_mi
+ max_mi = [0]*self.dim
+ for i in range(self.dim):
+ max_mi[i] = max(mi[i] for mi in
+ src_expansion.get_coefficient_identifiers())
+ max_mi[i] += max(mi[i] for mi in
+ self.get_coefficient_identifiers())
+
# Create a expansion terms wrangler for derivatives up to order
# (tgt order)+(src order).
tgtplusderiv_exp_terms_wrangler = \
src_expansion.expansion_terms_wrangler.copy(
- order=self.order + src_expansion.order)
+ order=self.order + src_expansion.order, max_mi=tuple(max_mi))
tgtplusderiv_coeff_ids = \
tgtplusderiv_exp_terms_wrangler.get_coefficient_identifiers()
tgtplusderiv_full_coeff_ids = \
diff --git a/sumpy/tools.py b/sumpy/tools.py
index 91c7605060bd46d071093a9c6131f7a71a7f89d7..8916e33c358702628d97cb03889ecb0f7aa861ec 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -401,48 +401,55 @@ def my_syntactic_subs(expr, subst_dict):
return expr
-# Source: http://scipy-cookbook.readthedocs.io/items/RankNullspace.html
-# Author: SciPy Developers
-# License: BSD-3-Clause
-
-def nullspace(A, atol=1e-13, rtol=0): # noqa:N803
- """Compute an approximate basis for the nullspace of A.
-
- The algorithm used by this function is based on the singular value
- decomposition of `A`.
-
- Parameters
- ----------
- A : ndarray
- A should be at most 2-D. A 1-D array with length k will be treated
- as a 2-D with shape (1, k)
- atol : float
- The absolute tolerance for a zero singular value. Singular values
- smaller than `atol` are considered to be zero.
- rtol : float
- The relative tolerance. Singular values less than rtol*smax are
- considered to be zero, where smax is the largest singular value.
-
- If both `atol` and `rtol` are positive, the combined tolerance is the
- maximum of the two; that is::
- tol = max(atol, rtol * smax)
- Singular values smaller than `tol` are considered to be zero.
-
- Return value
- ------------
- ns : ndarray
- If `A` is an array with shape (m, k), then `ns` will be an array
- with shape (k, n), where n is the estimated dimension of the
- nullspace of `A`. The columns of `ns` are a basis for the
- nullspace; each element in numpy.dot(A, ns) will be approximately
- zero.
- """
- from numpy.linalg import svd
- A = np.atleast_2d(A) # noqa:N806
- u, s, vh = svd(A)
- tol = max(atol, rtol * s[0])
- nnz = (s >= tol).sum()
- ns = vh[nnz:].conj().T
- return ns
+def rref(mat):
+ rows = len(mat)
+ cols = len(mat[0])
+ col = 0
+ pivot_cols = []
+
+ for row in range(rows):
+ if col >= cols:
+ break
+ i = row
+ while mat[i][col] == 0:
+ i += 1
+ if i == rows:
+ i = row
+ col += 1
+ if col == cols:
+ return mat, pivot_cols
+
+ pivot_cols.append(col)
+ mat[i], mat[row] = mat[row], mat[i]
+
+ piv = mat[row][col]
+ for c in range(col, cols):
+ mat[row][c] /= piv
+
+ for r in range(rows):
+ if r == row:
+ continue
+ piv = mat[r][col]
+ for c in range(col, cols):
+ mat[r][c] -= piv * mat[row][c]
+ col += 1
+ return mat, pivot_cols
+
+def nullspace(m):
+ m2 = [[sym.sympify(col) for col in row] for row in m]
+ mat, pivot_cols = rref(m2)
+ cols = len(mat[0])
+
+ free_vars = [i for i in range(cols) if i not in pivot_cols]
+
+ n = []
+ for free_var in free_vars:
+ vec = [0]*cols
+ vec[free_var] = 1
+ for piv_row, piv_col in enumerate(pivot_cols):
+ for pos in pivot_cols[piv_row+1:] + [free_var]:
+ vec[piv_col] -= mat[piv_row][pos]
+ n.append(vec)
+ return sym.Matrix(n).T
# vim: fdm=marker