Add orthonormal requirement

firedrakeproject · Aug 28, 2024 · 21f1bf1 · 21f1bf1
1 parent dfd1347
commit 21f1bf1
Show file tree

Hide file tree

Showing 4 changed files with 12 additions and 14 deletions.
diff --git a/FIAT/expansions.py b/FIAT/expansions.py
@@ -327,7 +327,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None):
         sd = self.ref_el.get_spatial_dimension()
 
         # Always return 1 for n=0 to make regression tests pass
-        scale = 1.0 if n == 0 and len(self.affine_mappings) == 1 else self.get_scale(cell=cell)
+        scale = 1.0 if n == 0 and (self.scale is None) else self.get_scale(cell=cell)
         phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv,
                                  scale, variant=self.variant)
         if self.continuity == "C0":

diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py
@@ -22,7 +22,7 @@ def NedelecSpace2D(ref_el, degree):
         raise Exception("NedelecSpace2D requires 2d reference element")
 
     k = degree - 1
-    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,))
+    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal")
 
     dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1)
     dimPk = expansions.polynomial_dimension(ref_el, k)
@@ -32,7 +32,7 @@ def NedelecSpace2D(ref_el, degree):
                                   for i in range(sd))))
     vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices)
 
-    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1)
+    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal")
     PkH = Pkp1.take(list(range(dimPkm1, dimPk)))
 
     Q = create_quadrature(ref_el, 2 * (k + 1))
@@ -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
@@ -169,10 +169,10 @@ 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.array(list(A.coeffs) + list(B.coeffs))
     assert A.get_reference_element() == B.get_reference_element()
     new_coeffs = construct_new_coeffs(A.get_reference_element(), A, B)
     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])
@@ -182,7 +182,9 @@ def polynomial_set_union_normalized(A, B):
         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)
@@ -206,13 +208,11 @@ def construct_new_coeffs(ref_el, A, B):
         higher = A if A.degree > B.degree else B
         lower = B if A.degree > B.degree else A
 
-        # dimHigher = expansions.polynomial_dimension(ref_el, higher.degree)
-        # dimLower = expansions.polynomial_dimension(ref_el, lower.degree)
-
         try:
             sd = lower.get_shape()[0]
         except IndexError:
             sd = 1
+
         embedded_coeffs = []
         diff = higher.coeffs.shape[-1] - lower.coeffs.shape[-1]
         for coeff in lower.coeffs:
@@ -224,11 +224,9 @@ def construct_new_coeffs(ref_el, A, B):
             else:
                 embedded_coeffs.append(numpy.append(coeff, [0 for i in range(diff)]))
         embedded_coeffs = numpy.array(embedded_coeffs)
-
         new_coeffs = numpy.array(list(embedded_coeffs) + list(higher.coeffs))
     else:
         new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs))
-
     return new_coeffs
 
 

diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py
@@ -20,7 +20,7 @@ def RTSpace(ref_el, degree):
     sd = ref_el.get_spatial_dimension()
 
     k = degree - 1
-    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,))
+    vec_Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, (sd,), scale="orthonormal")
 
     dimPkp1 = expansions.polynomial_dimension(ref_el, k + 1)
     dimPk = expansions.polynomial_dimension(ref_el, k)
@@ -30,7 +30,7 @@ def RTSpace(ref_el, degree):
                                   for i in range(sd))))
     vec_Pk_from_Pkp1 = vec_Pkp1.take(vec_Pk_indices)
 
-    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1)
+    Pkp1 = polynomial_set.ONPolynomialSet(ref_el, k + 1, scale="orthonormal")
     PkH = Pkp1.take(list(range(dimPkm1, dimPk)))
 
     Q = create_quadrature(ref_el, 2 * (k + 1))
@@ -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,