Merge branch 'master' into pbrubeck/gen-quad

FIAT/__init__.py

@@ -13,7 +13,7 @@
 from FIAT.hct import HsiehCloughTocher
 from FIAT.alfeld_sorokina import AlfeldSorokina
 from FIAT.arnold_qin import ArnoldQin
-from FIAT.guzman_neilan import GuzmanNeilan
+from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanH1div
 from FIAT.christiansen_hu import ChristiansenHu
 from FIAT.johnson_mercier import JohnsonMercier
 from FIAT.brezzi_douglas_marini import BrezziDouglasMarini
@@ -88,7 +88,9 @@
  "Alfeld-Sorokina": AlfeldSorokina,
  "Arnold-Qin": ArnoldQin,
  "Christiansen-Hu": ChristiansenHu,
- "Guzman-Neilan": GuzmanNeilan,
+ "Guzman-Neilan 1st kind H1": GuzmanNeilanFirstKindH1,
+ "Guzman-Neilan 2nd kind H1": GuzmanNeilanSecondKindH1,
+ "Guzman-Neilan H1(div)": GuzmanNeilanH1div,
  "Johnson-Mercier": JohnsonMercier,
  "Lagrange": Lagrange,
  "Kong-Mulder-Veldhuizen": KongMulderVeldhuizen,

FIAT/alfeld_sorokina.py

@@ -7,18 +7,16 @@
 # Written by Pablo D. Brubeck ([email protected]), 2024
 
 from FIAT import finite_element, dual_set, polynomial_set
-from FIAT.functional import ComponentPointEvaluation, PointDivergence, IntegralMoment
+from FIAT.functional import ComponentPointEvaluation, PointDivergence
 from FIAT.quadrature_schemes import create_quadrature
-from FIAT.quadrature import FacetQuadratureRule
 from FIAT.macro import CkPolynomialSet, AlfeldSplit
 
 import numpy
 
 
 def AlfeldSorokinaSpace(ref_el, degree):
- """Return a vector-valued C^0 PolynomialSet on an Alfeld split whose
- divergence is also C^0. This works on any simplex and for all
- polynomial degrees."""
+ """Return a vector-valued C0 PolynomialSet on an Alfeld split with C0
+ divergence. This works on any simplex and for all polynomial degrees."""
  ref_complex = AlfeldSplit(ref_el)
  sd = ref_complex.get_spatial_dimension()
  C0 = CkPolynomialSet(ref_complex, degree, order=0, shape=(sd,), variant="bubble")
@@ -51,53 +49,35 @@ def AlfeldSorokinaSpace(ref_el, degree):
 
 
 class AlfeldSorokinaDualSet(dual_set.DualSet):
- def __init__(self, ref_el, degree, variant=None):
+ def __init__(self, ref_el, degree):
  if degree != 2:
- raise NotImplementedError("Alfeld-Sorokina only defined for degree = 2")
- if variant is None:
- variant = "integral"
- if variant not in {"integral", "point"}:
- raise ValueError(f"Invalid variant {variant}")
+ raise NotImplementedError(f"{type(self).__name__} only defined for degree = 2")
 
  top = ref_el.get_topology()
  sd = ref_el.get_spatial_dimension()
  entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
 
  nodes = []
- dims = (0, 1) if variant == "point" else (0,)
- for dim in dims:
+ for dim in sorted(top):
  for entity in sorted(top[dim]):
  cur = len(nodes)
+
+ dpts = ref_el.make_points(dim, entity, degree-1)
+ nodes.extend(PointDivergence(ref_el, pt) for pt in dpts)
+
  pts = ref_el.make_points(dim, entity, degree)
- if dim == 0:
- pt, = pts
- nodes.append(PointDivergence(ref_el, pt))
  nodes.extend(ComponentPointEvaluation(ref_el, k, (sd,), pt)
  for pt in pts for k in range(sd))
  entity_ids[dim][entity].extend(range(cur, len(nodes)))
 
- if variant == "integral":
- # Edge moments against quadratic Legendre (mean-free bubble)
- dim = 1
- facet = ref_el.construct_subelement(dim)
- Q = create_quadrature(facet, degree+dim+1)
- f_at_qpts = facet.compute_bubble(Q.get_points())
- f_at_qpts -= numpy.dot(f_at_qpts, Q.get_weights()) / facet.volume()
- for entity in sorted(top[dim]):
- cur = len(nodes)
- Q_mapped = FacetQuadratureRule(ref_el, dim, entity, Q)
- detJ = Q_mapped.jacobian_determinant()
- phi = f_at_qpts / detJ
- nodes.extend(IntegralMoment(ref_el, Q_mapped, phi, comp=k, shp=(sd,))
- for k in range(sd))
- entity_ids[dim][entity].extend(range(cur, len(nodes)))
-
  super().__init__(nodes, ref_el, entity_ids)
 
 
 class AlfeldSorokina(finite_element.CiarletElement):
- """The Alfeld-Sorokina C^0 quadratic macroelement with C^0 divergence. This element
- belongs to a Stokes complex, and is paired with Lagrange(ref_el, 1, variant="alfeld")."""
+ """The Alfeld-Sorokina C0 quadratic macroelement with C0 divergence.
+
+ This element belongs to a Stokes complex, and is paired with CG1(Alfeld).
+ """
  def __init__(self, ref_el, degree=2):
  dual = AlfeldSorokinaDualSet(ref_el, degree)
  poly_set = AlfeldSorokinaSpace(ref_el, degree)

FIAT/arnold_qin.py

@@ -16,8 +16,8 @@
 
 def ArnoldQinSpace(ref_el, degree, reduced=False):
  """Return a basis for the Arnold-Qin space
- curl(HCT-red) + P_0 x if reduced = True, and
- curl(HCT) + P_0 x if reduced = False."""
+ curl(HCT-red) + P0 x if reduced = True, and
+ curl(HCT) + P0 x if reduced = False."""
  if ref_el.get_shape() != TRIANGLE:
  raise ValueError("Arnold-Qin only defined on triangles")
  if degree != 2:
@@ -56,16 +56,16 @@ def ArnoldQinSpace(ref_el, degree, reduced=False):
 
 
 class ArnoldQin(finite_element.CiarletElement):
- """The Arnold-Qin C^0(Alfeld) quadratic macroelement with divergence in P0.
+ """The Arnold-Qin C0(Alfeld) quadratic macroelement with divergence in P0.
  This element belongs to a Stokes complex, and is paired with unsplit DG0."""
  def __init__(self, ref_el, degree=2, reduced=False):
  poly_set = ArnoldQinSpace(ref_el, degree)
  if reduced:
- subdegree = 1
+ order = 1
  mapping = "contravariant piola"
  else:
- subdegree = degree
+ order = degree
  mapping = "affine"
- dual = BernardiRaugelDualSet(ref_el, degree, subdegree)
+ dual = BernardiRaugelDualSet(ref_el, order, degree=degree)
  formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form
  super().__init__(poly_set, dual, degree, formdegree, mapping=mapping)

FIAT/bernardi_raugel.py

@@ -11,92 +11,98 @@
 
 from FIAT import finite_element, dual_set, polynomial_set, expansions
 from FIAT.functional import ComponentPointEvaluation, FrobeniusIntegralMoment
-from FIAT.quadrature_schemes import create_quadrature
+from FIAT.hierarchical import make_dual_bubbles
 from FIAT.quadrature import FacetQuadratureRule
+
 import numpy
+import math
 
 
-def ExtendedBernardiRaugelSpace(ref_el, subdegree):
- r"""Return a basis for the extended Bernardi-Raugel space.
- P_k^d + (P_{dim} \ P_{dim-1})^d"""
+def BernardiRaugelSpace(ref_el, order):
+ """Return a basis for the extended Bernardi-Raugel space: (Pk + FacetBubble)^d."""
  sd = ref_el.get_spatial_dimension()
- if subdegree >= sd:
- raise ValueError("The Bernardi-Raugel space is only defined for subdegree < dim")
+ if order > sd:
+ raise ValueError("The Bernardi-Raugel space is only defined for order <= dim")
  Pd = polynomial_set.ONPolynomialSet(ref_el, sd, shape=(sd,), scale=1, variant="bubble")
  dimPd = expansions.polynomial_dimension(ref_el, sd, continuity="C0")
  entity_ids = expansions.polynomial_entity_ids(ref_el, sd, continuity="C0")
+
+ slices = {dim: slice(math.comb(order-1, dim)) for dim in range(order)}
+ slices[sd-1] = slice(None)
  ids = [i + j * dimPd
- for dim in (*tuple(range(subdegree)), sd-1)
+ for dim in slices
  for f in sorted(entity_ids[dim])
- for i in entity_ids[dim][f]
+ for i in entity_ids[dim][f][slices[dim]]
  for j in range(sd)]
  return Pd.take(ids)
 
 
 class BernardiRaugelDualSet(dual_set.DualSet):
  """The Bernardi-Raugel dual set."""
- def __init__(self, ref_complex, degree, subdegree=1, reduced=False):
- ref_el = ref_complex.get_parent() or ref_complex
+ def __init__(self, ref_el, order=1, degree=None, reduced=False, ref_complex=None):
+ if ref_complex is None:
+ ref_complex = ref_el
  sd = ref_el.get_spatial_dimension()
- if subdegree > sd:
- raise ValueError("The Bernardi-Raugel dual is only defined for subdegree <= dim")
+ if degree is None:
+ degree = sd
+ if order > sd:
+ raise ValueError(f"{type(self).__name__} is only defined for order <= dim")
  top = ref_el.get_topology()
  entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)}
 
- # Point evaluation at lattice points
  nodes = []
- for dim in range(subdegree):
- for entity in sorted(top[dim]):
- cur = len(nodes)
- pts = ref_el.make_points(dim, entity, subdegree)
- nodes.extend(ComponentPointEvaluation(ref_el, comp, (sd,), pt)
- for pt in pts for comp in range(sd))
- entity_ids[dim][entity].extend(range(cur, len(nodes)))
-
- if subdegree < sd:
- # Face moments of normal/tangential components against mean-free bubbles
- facet = ref_complex.construct_subcomplex(sd-1)
- Q = create_quadrature(facet, 2*degree)
- if degree == 1 and facet.is_macrocell():
- P = polynomial_set.ONPolynomialSet(facet, degree, scale=1, variant="bubble")
- f_at_qpts = P.tabulate(Q.get_points())[(0,)*(sd-1)][-1]
- else:
- ref_facet = facet.get_parent() or facet
- f_at_qpts = ref_facet.compute_bubble(Q.get_points())
- f_at_qpts -= numpy.dot(f_at_qpts, Q.get_weights()) / facet.volume()
-
- Qs = {f: FacetQuadratureRule(ref_el, sd-1, f, Q)
- for f in sorted(top[sd-1])}
-
- thats = {f: ref_el.compute_tangents(sd-1, f)
- for f in sorted(top[sd-1])}
+ if order > 0:
+ # Point evaluation at lattice points
+ for dim in sorted(top):
+ for entity in sorted(top[dim]):
+ cur = len(nodes)
+ pts = ref_el.make_points(dim, entity, order)
+ nodes.extend(ComponentPointEvaluation(ref_el, comp, (sd,), pt)
+ for pt in pts for comp in range(sd))
+ entity_ids[dim][entity].extend(range(cur, len(nodes)))
+
+ if order < sd:
+ # Face moments of normal/tangential components against dual bubbles
+ ref_facet = ref_complex.construct_subcomplex(sd-1)
+ codim = sd-1 if degree == 1 and ref_facet.is_macrocell() else 0
+ Q, phis = make_dual_bubbles(ref_facet, degree, codim=codim)
+ f_at_qpts = phis[-1]
+ if codim != 0:
+ f_at_qpts -= numpy.dot(f_at_qpts, Q.get_weights()) / ref_facet.volume()
+
+ interior_facets = ref_el.get_interior_facets(sd-1) or ()
+ facets = list(set(top[sd-1]) - set(interior_facets))
+
+ Qs = {f: FacetQuadratureRule(ref_el, sd-1, f, Q) for f in facets}
+ thats = {f: ref_el.compute_tangents(sd-1, f) for f in facets}
 
  R = numpy.array([[0, 1], [-1, 0]])
  ndir = 1 if reduced else sd
  for i in range(ndir):
- for f, Q_mapped in Qs.items():
+ for f in sorted(facets):
  cur = len(nodes)
  if i == 0:
  udir = numpy.dot(R, *thats[f]) if sd == 2 else numpy.cross(*thats[f])
  else:
  udir = thats[f][i-1]
- detJ = Q_mapped.jacobian_determinant()
+ detJ = Qs[f].jacobian_determinant()
  phi_at_qpts = udir[:, None] * f_at_qpts[None, :] / detJ
- nodes.append(FrobeniusIntegralMoment(ref_el, Q_mapped, phi_at_qpts))
+ nodes.append(FrobeniusIntegralMoment(ref_el, Qs[f], phi_at_qpts))
  entity_ids[sd-1][f].extend(range(cur, len(nodes)))
-
  super().__init__(nodes, ref_el, entity_ids)
 
 
 class BernardiRaugel(finite_element.CiarletElement):
- """The Bernardi-Raugel extended element."""
- def __init__(self, ref_el, degree=None, subdegree=1):
- sd = ref_el.get_spatial_dimension()
- if degree is None:
- degree = sd
- if degree != sd:
- raise ValueError("Bernardi-Raugel only defined for degree = dim")
- poly_set = ExtendedBernardiRaugelSpace(ref_el, subdegree)
- dual = BernardiRaugelDualSet(ref_el, degree, subdegree=subdegree)
- formdegree = sd - 1 # (n-1)-form
+ """The Bernardi-Raugel (extended) element.
+
+ This element does not belong to a Stokes complex, but can be paired with
+ DG_{k-1}. This pair is inf-sup stable, but only weakly divergence-free.
+ """
+ def __init__(self, ref_el, order=1):
+ degree = ref_el.get_spatial_dimension()
+ if order >= degree:
+ raise ValueError(f"{type(self).__name__} only defined for order < dim")
+ poly_set = BernardiRaugelSpace(ref_el, order)
+ dual = BernardiRaugelDualSet(ref_el, order, degree=degree)
+ formdegree = 0
  super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola")

FIAT/christiansen_hu.py

@@ -40,7 +40,7 @@ def ChristiansenHuSpace(ref_el, degree, reduced=False):
 
  if not reduced:
  # Compute the primal basis via Vandermonde and extract the facet bubbles
- dual = BernardiRaugelDualSet(ref_complex, degree, reduced=True)
+ dual = BernardiRaugelDualSet(ref_el, degree, degree=degree, ref_complex=ref_complex, reduced=True)
  dualmat = dual.to_riesz(C0)
  V = numpy.tensordot(dualmat, coeffs, axes=((1, 2), (1, 2)))
  coeffs = numpy.tensordot(numpy.linalg.inv(V.T), coeffs, axes=(-1, 0))
@@ -73,6 +73,6 @@ def __init__(self, ref_el, degree=1):
  raise ValueError("Christiansen-Hu only defined for degree = 1")
  poly_set = ChristiansenHuSpace(ref_el, degree)
  ref_complex = poly_set.get_reference_element()
- dual = BernardiRaugelDualSet(ref_complex, degree)
+ dual = BernardiRaugelDualSet(ref_el, degree, degree=degree, ref_complex=ref_complex)
  formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form
  super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola")

FIAT/dual_set.py

@@ -253,7 +253,7 @@ def make_entity_closure_ids(ref_el, entity_ids):
 
 def lexsort_nodes(ref_el, nodes, entity=None, offset=0):
  """Sort PointEvaluation nodes in lexicographical ordering."""
- if len(nodes) > 1 and all(isinstance(node, functional.PointEvaluation) for node in nodes):
+ if len(nodes) > 1:
  pts = []
  for node in nodes:
  pt, = node.get_point_dict()
@@ -270,17 +270,29 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations):
  parent_cell = ref_el.get_parent()
  if parent_cell is None:
  return nodes, ref_el, entity_ids, entity_permutations
- parent_nodes = []
  parent_ids = {}
  parent_permutations = None
-
  parent_to_children = ref_el.get_parent_to_children()
- for dim in sorted(parent_to_children):
- parent_ids[dim] = {}
- for entity in sorted(parent_to_children[dim]):
- cur = len(parent_nodes)
- for child_dim, child_entity in parent_to_children[dim][entity]:
- parent_nodes.extend(nodes[i] for i in entity_ids[child_dim][child_entity])
- ids = lexsort_nodes(parent_cell, parent_nodes[cur:], entity=(dim, entity), offset=cur)
- parent_ids[dim][entity] = ids
+
+ if all(isinstance(node, functional.PointEvaluation) for node in nodes):
+ # Merge Lagrange dual with lexicographical reordering
+ parent_nodes = []
+ for dim in sorted(parent_to_children):
+ parent_ids[dim] = {}
+ for entity in sorted(parent_to_children[dim]):
+ cur = len(parent_nodes)
+ for child_dim, child_entity in parent_to_children[dim][entity]:
+ parent_nodes.extend(nodes[i] for i in entity_ids[child_dim][child_entity])
+ ids = lexsort_nodes(parent_cell, parent_nodes[cur:], entity=(dim, entity), offset=cur)
+ parent_ids[dim][entity] = ids
+ else:
+ # Merge everything else with the same node ordering
+ parent_nodes = nodes
+ for dim in sorted(parent_to_children):
+ parent_ids[dim] = {}
+ for entity in sorted(parent_to_children[dim]):
+ parent_ids[dim][entity] = []
+ for child_dim, child_entity in parent_to_children[dim][entity]:
+ parent_ids[dim][entity].extend(entity_ids[child_dim][child_entity])
+
  return parent_nodes, parent_cell, parent_ids, parent_permutations

FIAT/expansions.py

@@ -512,6 +512,11 @@ def tabulate_jet(self, n, pts, order=1):
  data.append(vr.transpose((r, r+1) + tuple(range(r))))
  return data
 
+ def __eq__(self, other):
+ return (type(self) is type(other) and
+ self.ref_el == other.ref_el and
+ self.continuity == other.continuity)
+
 
 class PointExpansionSet(ExpansionSet):
  """Evaluates the point basis on a point reference element."""

FIAT/finite_element.py

@@ -132,7 +132,7 @@ class CiarletElement(FiniteElement):
  def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref_complex=None):
  ref_el = dual.get_reference_element()
  ref_complex = ref_complex or poly_set.get_reference_element()
- super(CiarletElement, self).__init__(ref_el, dual, order, formdegree, mapping, ref_complex)
+ super().__init__(ref_el, dual, order, formdegree, mapping, ref_complex)
 
  # build generalized Vandermonde matrix
  old_coeffs = poly_set.get_coeffs()