Extend PolySet coefficients (#130)

* Add extension of coefficients to allow union of polysets with different degrees

* remove orthonormal as seems to be no longer necessary

* Remove further orthonormals

Author: indiamai

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,

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.

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,

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)