From 98fc7595b5d4ef048b40f04db35b245979fa5b34 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 6 Nov 2024 23:30:29 +0000 Subject: [PATCH] Fix serendipitous supersmoothness assumption (#102) * Fix serendipitous supersmoothness assumption * unify svd computations --- FIAT/alfeld_sorokina.py | 6 ++---- FIAT/christiansen_hu.py | 6 ++---- FIAT/macro.py | 25 ++++++++++--------------- FIAT/polynomial_set.py | 29 +++++++++++++---------------- 4 files changed, 27 insertions(+), 39 deletions(-) diff --git a/FIAT/alfeld_sorokina.py b/FIAT/alfeld_sorokina.py index 5d3244b78..e09f26d17 100644 --- a/FIAT/alfeld_sorokina.py +++ b/FIAT/alfeld_sorokina.py @@ -41,10 +41,8 @@ def AlfeldSorokinaSpace(ref_el, degree): if len(rows) > 0: dual_mat = numpy.vstack(rows) - _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) - tol = sig[0] * 1E-10 - num_sv = len([s for s in sig if abs(s) > tol]) - coeffs = numpy.tensordot(vt[num_sv:], coeffs, axes=(-1, 0)) + nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) + coeffs = numpy.tensordot(nsp, coeffs, axes=(-1, 0)) return polynomial_set.PolynomialSet(ref_complex, degree, degree, expansion_set, coeffs) diff --git a/FIAT/christiansen_hu.py b/FIAT/christiansen_hu.py index 476ea914d..64fe89659 100644 --- a/FIAT/christiansen_hu.py +++ b/FIAT/christiansen_hu.py @@ -28,10 +28,8 @@ def ChristiansenHuSpace(ref_el, degree, reduced=False): tab = C0.tabulate(Q.get_points(), 1) divC0 = sum(tab[alpha][:, alpha.index(1), :] for alpha in tab if sum(alpha) == 1) - _, sig, vt = numpy.linalg.svd(divC0.T, full_matrices=True) - tol = sig[0] * 1E-10 - num_sv = len([s for s in sig if abs(s) > tol]) - coeffs = numpy.tensordot(vt[num_sv:], C0.get_coeffs(), axes=(-1, 0)) + nsp = polynomial_set.spanning_basis(divC0.T, nullspace=True) + coeffs = numpy.tensordot(nsp, C0.get_coeffs(), axes=(-1, 0)) verts = numpy.array(ref_complex.get_vertices()) WT = verts[-1] diff --git a/FIAT/macro.py b/FIAT/macro.py index f48298c70..77ab89382 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -429,25 +429,21 @@ def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs): if len(rows) > 0: dual_mat = numpy.vstack(rows) - _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) - tol = sig[0] * 1E-10 - num_sv = len([s for s in sig if abs(s) > tol]) - coeffs = vt[num_sv:] + coeffs = polynomial_set.spanning_basis(dual_mat, nullspace=True) else: coeffs = numpy.eye(expansion_set.get_num_members(degree)) - if vorder > order + 1: - # Impose C^vorder super-smoothness at interior vertices - # C^order automatically gives C^{order+1} at the interior vertex + # Impose C^vorder super-smoothness at interior vertices + # C^order automatically gives C^{order+dim-1} at the interior vertex + sorder = order + sd - 1 + if vorder > sorder: verts = ref_el.get_vertices() points = [verts[i] for i in ref_el.get_interior_facets(0)] jumps = expansion_set.tabulate_jumps(degree, points, order=vorder) - for r in range(order+2, vorder+1): + for r in range(sorder+1, vorder+1): dual_mat = numpy.dot(numpy.vstack(jumps[r].T), coeffs.T) - _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) - tol = sig[0] * 1E-10 - num_sv = len([s for s in sig if abs(s) > tol]) - coeffs = numpy.dot(vt[num_sv:], coeffs) + nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) + coeffs = numpy.dot(nsp, coeffs) if shape != (): m, n = coeffs.shape @@ -498,8 +494,7 @@ def __init__(self, ref_el, degree, order=0, **kwargs): if len(rows) > 0: dual_mat = numpy.vstack(rows) - _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) - num_sv = len([s for s in sig if abs(s) > 1.e-10]) - coeffs = numpy.tensordot(vt[num_sv:], coeffs, axes=(1, 0)) + nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) + coeffs = numpy.tensordot(nsp, coeffs, axes=(1, 0)) super().__init__(ref_el, degree, degree, expansion_set, coeffs) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 579ec3614..3a3e9ce3a 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -162,28 +162,25 @@ def form_matrix_product(mats, alpha): return result +def spanning_basis(A, nullspace=False, rtol=1e-10): + """Construct a basis that spans the rows of A via SVD. + """ + Aflat = A.reshape(A.shape[0], -1) + u, sig, vt = numpy.linalg.svd(Aflat, full_matrices=True) + atol = rtol * (sig[0] + 1) + num_sv = len([s for s in sig if abs(s) > atol]) + basis = vt[num_sv:] if nullspace else vt[:num_sv] + return numpy.reshape(basis, (-1, *A.shape[1:])) + + def polynomial_set_union_normalized(A, B): """Given polynomial sets A and B, constructs a new polynomial set whose span is the same as that of span(A) union span(B). It may not contain any of the same members of the set, as we construct a span via SVD. """ - new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) - func_shape = new_coeffs.shape[1:] - if len(func_shape) == 1: - (u, sig, vt) = numpy.linalg.svd(new_coeffs) - num_sv = len([s for s in sig if abs(s) > 1.e-10]) - coeffs = vt[:num_sv] - else: - new_shape0 = new_coeffs.shape[0] - new_shape1 = numpy.prod(func_shape) - newshape = (new_shape0, new_shape1) - nc = numpy.reshape(new_coeffs, newshape) - (u, sig, vt) = numpy.linalg.svd(nc, 1) - num_sv = len([s for s in sig if abs(s) > 1.e-10]) - - coeffs = numpy.reshape(vt[:num_sv], (num_sv,) + func_shape) - + new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) + coeffs = spanning_basis(new_coeffs) return PolynomialSet(A.get_reference_element(), A.get_degree(), A.get_embedded_degree(),