fix JM dof ordering

firedrakeproject · Oct 27, 2024 · a4cdfec · a4cdfec
1 parent 2075f9b
commit a4cdfec
Show file tree

Hide file tree

Showing 3 changed files with 5 additions and 5 deletions.
diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py
@@ -9,7 +9,7 @@
 # SPDX-License-Identifier:    LGPL-3.0-or-later
 
 
-from FIAT import finite_element, dual_set, expansions, polynomial_set
+from FIAT import finite_element, dual_set, polynomial_set
 from FIAT.reference_element import TRIANGLE
 from FIAT.quadrature_schemes import create_quadrature
 from FIAT.functional import (ComponentPointEvaluation,
@@ -109,9 +109,9 @@ def __init__(self, ref_el, degree=3):
                      for phi in phis for i in range(sd) for j in range(i, sd))
 
         # constraint dofs: moments of divergence against P_{k-1} \ P_{k-2}
-        dimPkm1 = expansions.polynomial_dimension(ref_el, degree-1)
-        dimPkm2 = expansions.polynomial_dimension(ref_el, degree-2)
         P = polynomial_set.ONPolynomialSet(ref_el, degree-1, shape=(sd,))
+        dimPkm1 = P.expansion_set.get_num_members(degree-1)
+        dimPkm2 = P.expansion_set.get_num_members(degree-2)
         PH = P.take([i + j * dimPkm1 for j in range(sd) for i in range(dimPkm2, dimPkm1)])
         phis = PH.tabulate(Q.get_points())[(0,)*sd]
         nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis)

diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py
@@ -84,4 +84,4 @@ def __init__(self, ref_el, degree=3, variant=None):
         poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree)
         dual = HuZhangDual(ref_el, degree, variant, qdegree)
         formdegree = ref_el.get_spatial_dimension() - 1
-        super().__init__(poly_set, dual, degree, formdegree=formdegree, mapping="double contravariant piola")
+        super().__init__(poly_set, dual, degree, formdegree, mapping="double contravariant piola")
diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py
@@ -44,7 +44,7 @@ def __init__(self, ref_complex, degree, variant=None):
             P = polynomial_set.ONPolynomialSet(ref_el, degree-1)
             phis = P.tabulate(Q.get_points())[(0,) * sd]
             nodes.extend(IntegralMoment(ref_el, Q, phi, comp=(i, j))
-                         for j in range(sd) for i in range(j+1) for phi in phis)
+                         for phi in phis for i in range(sd) for j in range(i, sd))
 
         entity_ids[sd][0].extend(range(cur, len(nodes)))