diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index c92fa06be..6cf122772 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -43,8 +43,8 @@ def NedelecSpace2D(ref_el, degree): CrossX = numpy.dot(numpy.array([[0.0, 1.0], [-1.0, 0.0]]), Qpts.T) PkHCrossX_at_Qpts = PkH_at_Qpts[:, None, :] * CrossX[None, :, :] - PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), Pkp1_at_Qpts.T) - + PkHCrossX_coeffs = numpy.dot(numpy.multiply(PkHCrossX_at_Qpts, Qwts), + Pkp1_at_Qpts.T) PkHcrossX = polynomial_set.PolynomialSet(ref_el, k + 1, k + 1, diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 08e62b8d8..9686cf245 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -176,15 +176,49 @@ def polynomial_set_union_normalized(A, B): not contain any of the same members of the set, as we construct a span via SVD. """ - new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) + assert A.get_reference_element() == B.get_reference_element() + new_coeffs = construct_new_coeffs(A.get_reference_element(), A, B) + + deg = max(A.get_degree(), B.get_degree()) + em_deg = max(A.get_embedded_degree(), B.get_embedded_degree()) coeffs = spanning_basis(new_coeffs) return PolynomialSet(A.get_reference_element(), - A.get_degree(), - A.get_embedded_degree(), + deg, + em_deg, A.get_expansion_set(), coeffs) +def construct_new_coeffs(ref_el, A, B): + """ + Constructs new coefficients for the set union of A and B + If A and B are discontinuous and do not have the same degree the smaller one + is extended to match the larger. + + This does not handle the case that A and B have continuity but not the same degree. + """ + + if A.get_expansion_set().continuity != B.get_expansion_set().continuity: + raise ValueError("Continuity of expansion sets does not match.") + + if A.get_embedded_degree() != B.get_embedded_degree() and A.get_expansion_set().continuity is None: + higher = A if A.get_embedded_degree() > B.get_embedded_degree() else B + lower = B if A.get_embedded_degree() > B.get_embedded_degree() else A + + diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1] + + # pad only the 0th axis with the difference in size + padding = [(0, 0) for i in range(len(lower.coeffs.shape) - 1)] + [(0, diff)] + embedded_coeffs = numpy.pad(lower.coeffs, padding) + + new_coeffs = numpy.concatenate((embedded_coeffs, higher.coeffs), axis=0) + elif A.get_embedded_degree() == B.get_embedded_degree(): + new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) + else: + raise NotImplementedError("Extending of coefficients is not implemented for PolynomialSets with continuity and different degrees") + return new_coeffs + + class ONSymTensorPolynomialSet(PolynomialSet): """Constructs an orthonormal basis for symmetric-tensor-valued polynomials on a reference element. diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index db0058d12..2c7d72c83 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -39,12 +39,12 @@ def RTSpace(ref_el, degree): # have to work on this through "tabulate" interface # first, tabulate PkH at quadrature points PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd] + Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd] x = Qpts.T PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :] PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T) - PkHx = polynomial_set.PolynomialSet(ref_el, k, k + 1, diff --git a/test/FIAT/unit/test_polynomial.py b/test/FIAT/unit/test_polynomial.py index 2982cf11e..c211fb815 100644 --- a/test/FIAT/unit/test_polynomial.py +++ b/test/FIAT/unit/test_polynomial.py @@ -21,6 +21,7 @@ from FIAT import expansions, polynomial_set, reference_element from FIAT.quadrature_schemes import create_quadrature +from itertools import chain @pytest.fixture(params=(1, 2, 3)) @@ -107,3 +108,45 @@ def test_bubble_duality(cell, degree): results = scale * numpy.dot(numpy.multiply(phi_dual, qwts), phi.T) assert numpy.allclose(results, numpy.diag(numpy.diag(results))) assert numpy.allclose(numpy.diag(results), 1.0) + + +@pytest.mark.parametrize("degree", [10]) +def test_union_of_polysets(cell, degree): + """ demonstrates that polysets don't need to have the same degree for union + using RT space as an example""" + + sd = cell.get_spatial_dimension() + k = degree + vecPk = polynomial_set.ONPolynomialSet(cell, degree, (sd,)) + + vec_Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, (sd,), scale="orthonormal") + + dimPkp1 = expansions.polynomial_dimension(cell, k + 1) + dimPk = expansions.polynomial_dimension(cell, k) + dimPkm1 = expansions.polynomial_dimension(cell, k - 1) + + vec_Pk_indices = list(chain(*(range(i * dimPkp1, i * dimPkp1 + dimPk) + for i in range(sd)))) + vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices) + + Pkp1 = polynomial_set.ONPolynomialSet(cell, k + 1, scale="orthonormal") + PkH = Pkp1.take(list(range(dimPkm1, dimPk))) + + Q = create_quadrature(cell, 2 * (k + 1)) + Qpts, Qwts = Q.get_points(), Q.get_weights() + + PkH_at_Qpts = PkH.tabulate(Qpts)[(0,) * sd] + Pkp1_at_Qpts = Pkp1.tabulate(Qpts)[(0,) * sd] + x = Qpts.T + PkHx_at_Qpts = PkH_at_Qpts[:, None, :] * x[None, :, :] + PkHx_coeffs = numpy.dot(numpy.multiply(PkHx_at_Qpts, Qwts), Pkp1_at_Qpts.T) + PkHx = polynomial_set.PolynomialSet(cell, k, k + 1, vec_Pkp1.get_expansion_set(), PkHx_coeffs) + + same_deg = polynomial_set.polynomial_set_union_normalized(vec_Pk_from_Pkp1, PkHx) + different_deg = polynomial_set.polynomial_set_union_normalized(vecPk, PkHx) + + Q = create_quadrature(cell, 2*(degree)) + Qpts, _ = Q.get_points(), Q.get_weights() + same_vals = same_deg.tabulate(Qpts)[(0,) * sd] + diff_vals = different_deg.tabulate(Qpts)[(0,) * sd] + assert numpy.allclose(same_vals - diff_vals, 0)