Skip to content

Harmonic Guzman-Neilan #232

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
pbrubeck wants to merge 2 commits into
base: main
Choose a base branch
from
Draft

Harmonic Guzman-Neilan #232

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9fa98bd
Harmonic Guzman-Neilan
pbrubeck Mar 16, 2026
12cf7d5
Merge branch 'main' into pbrubeck/gn-harmonic
pbrubeck Apr 28, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
48 changes: 36 additions & 12 deletions FIAT/guzman_neilan.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@
import math


def GuzmanNeilanSpace(ref_el, order, kind=1, reduced=False):
def GuzmanNeilanSpace(ref_el, order, kind=1, reduced=False, harmonic=False):
r"""Return a basis for the (extended) Guzman-Neilan H1 space.

Project the extended Bernardi-Raugel space (Pk + FacetBubble)^d
Expand All @@ -33,14 +33,15 @@ def GuzmanNeilanSpace(ref_el, order, kind=1, reduced=False):
:kwarg kind: kind = 1 gives Pk^d + GN bubbles,
kind = 2 gives C0 Pk(Alfeld)^d + GN bubbles.
:kwarg reduced: Include tangential bubbles if reduced = False.
:kwarg harmonic: Use Stokes-harmonic extension for the face bubbles.

:returns: a PolynomialSet basis for the Guzman-Neilan H1 space.
"""
sd = ref_el.get_spatial_dimension()
ref_complex = AlfeldSplit(ref_el)
C0 = polynomial_set.ONPolynomialSet(ref_complex, sd, shape=(sd,), scale=1, variant="bubble")
B = take_interior_bubbles(C0)
if sd > 2:
if sd > 2 and not harmonic:
B = modified_bubble_subspace(B)

K = ref_complex if kind == 2 else ref_el
Expand All @@ -59,12 +60,18 @@ def GuzmanNeilanSpace(ref_el, order, kind=1, reduced=False):

class GuzmanNeilanH1(finite_element.CiarletElement):
"""The Guzman-Neilan H1-conforming (extended) macroelement."""
def __init__(self, ref_el, order=1, kind=1, quad_scheme=None):
def __init__(self, ref_el, order=1, kind=1, variant=None, quad_scheme=None):
sd = ref_el.get_spatial_dimension()
if order >= sd:
raise ValueError(f"{type(self).__name__} is only defined for order < dim")
harmonic = False
if variant is not None:
if variant == "harmonic":
harmonic = True
else:
raise ValueError(f"{type(self).__name__} does not support variant {variant}")
degree = sd
poly_set = GuzmanNeilanSpace(ref_el, order, kind=kind)
poly_set = GuzmanNeilanSpace(ref_el, order, kind=kind, harmonic=harmonic)
ref_complex = poly_set.get_reference_element() if kind == 2 else ref_el
dual = BernardiRaugelDualSet(ref_complex, order, degree=degree, quad_scheme=quad_scheme)
formdegree = sd - 1 # (n-1)-form
Expand All @@ -80,8 +87,8 @@ class GuzmanNeilanFirstKindH1(GuzmanNeilanH1):

This element belongs to a Stokes complex, and is paired with unsplit DG_{k-1}.
"""
def __init__(self, ref_el, order=1, quad_scheme=None):
super().__init__(ref_el, order=order, kind=1, quad_scheme=quad_scheme)
def __init__(self, ref_el, order=1, variant=None, quad_scheme=None):
super().__init__(ref_el, order=order, kind=1, variant=variant, quad_scheme=quad_scheme)


class GuzmanNeilanSecondKindH1(GuzmanNeilanH1):
Expand All @@ -93,11 +100,11 @@ class GuzmanNeilanSecondKindH1(GuzmanNeilanH1):

This element belongs to a Stokes complex, and is paired with DG_{k-1}(Alfeld).
"""
def __init__(self, ref_el, order=1, quad_scheme=None):
super().__init__(ref_el, order=order, kind=2, quad_scheme=quad_scheme)
def __init__(self, ref_el, order=1, variant=None, quad_scheme=None):
super().__init__(ref_el, order=order, kind=2, variant=variant, quad_scheme=quad_scheme)


def GuzmanNeilanH1div(ref_el, degree=2, reduced=False, quad_scheme=None):
def GuzmanNeilanH1div(ref_el, degree=2, variant=None, reduced=False, quad_scheme=None):
"""The Guzman-Neilan H1(div)-conforming (extended) macroelement.

Reference element: a simplex of any dimension.
Expand All @@ -115,7 +122,7 @@ def GuzmanNeilanH1div(ref_el, degree=2, reduced=False, quad_scheme=None):
div_nodes = [i for i, node in enumerate(AS.dual_basis())
if len(node.deriv_dict) > 0]
AS = RestrictedElement(AS, indices=div_nodes)
GN = GuzmanNeilanH1(ref_el, order=order, quad_scheme=quad_scheme)
GN = GuzmanNeilanH1(ref_el, order=order, variant=variant, quad_scheme=quad_scheme)
return NodalEnrichedElement(AS, GN)


Expand All @@ -130,6 +137,11 @@ def div(U):
return sum(U[k][:, k.index(1), :] for k in U if sum(k) == 1)


def grad(U):
"""Compute the gradient from tabulation dict."""
return numpy.concatenate([U[k] for k in U if sum(k) == 1], axis=-2)


def take_interior_bubbles(P, degree=None):
"""Extract the interior bubbles up to the given degree from a complete PolynomialSet."""
ref_complex = P.get_reference_element()
Expand Down Expand Up @@ -202,10 +214,22 @@ def constant_div_projection(BR, C0, M, num_bubbles):

U = M.tabulate(qpts, 1)
X = BR.tabulate(qpts, 1)
# Invert the divergence
B = inner(P, div(U), qwts)
g = inner(P, div(X)[-num_bubbles:], qwts)
w = numpy.linalg.solve(B, g)
if B.shape[0] == B.shape[1]:
# Invert the divergence
w = numpy.linalg.solve(B, g)
else:
# Solve Stokes
n, m = B.shape
A = inner(grad(U), grad(U), qwts)
S = numpy.zeros((n + m, n + m))
S[:m, :m] = A
S[m:, :m] = B
S[:m, m:] = B.T
f = inner(grad(U), grad(X)[-num_bubbles:], qwts)
rhs = numpy.concatenate([f, g], axis=0)
w = numpy.linalg.solve(S, rhs)[:m]

# Add correction to BR bubbles
v = C0.tabulate(qpts)[(0,)*sd]
Expand Down
16 changes: 8 additions & 8 deletions finat/guzman_neilan.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,26 +6,26 @@

class GuzmanNeilanFirstKindH1(PiolaBubbleElement):
"""Pk^d enriched with Guzman-Neilan bubbles."""
def __init__(self, cell, order=1, quad_scheme=None):
def __init__(self, cell, order=1, variant=None, quad_scheme=None):
cite("GuzmanNeilan2018")
super().__init__(FIAT.GuzmanNeilanFirstKindH1(cell, order=order, quad_scheme=quad_scheme))
super().__init__(FIAT.GuzmanNeilanFirstKindH1(cell, order=order, variant=variant, quad_scheme=quad_scheme))


class GuzmanNeilanSecondKindH1(PiolaBubbleElement):
"""C0 Pk^d(Alfeld) enriched with Guzman-Neilan bubbles."""
def __init__(self, cell, order=1, quad_scheme=None):
def __init__(self, cell, order=1, variant=None, quad_scheme=None):
cite("GuzmanNeilan2018")
super().__init__(FIAT.GuzmanNeilanSecondKindH1(cell, order=order, quad_scheme=quad_scheme))
super().__init__(FIAT.GuzmanNeilanSecondKindH1(cell, order=order, variant=variant, quad_scheme=quad_scheme))


class GuzmanNeilanBubble(GuzmanNeilanFirstKindH1):
"""Modified Bernardi-Raugel bubbles that are C^0 P_dim(Alfeld) with constant divergence."""
def __init__(self, cell, degree=None, quad_scheme=None):
super().__init__(cell, order=0, quad_scheme=quad_scheme)
def __init__(self, cell, degree=None, variant=None, quad_scheme=None):
super().__init__(cell, order=0, variant=variant, quad_scheme=quad_scheme)


class GuzmanNeilanH1div(PiolaBubbleElement):
"""Alfeld-Sorokina nodally enriched with Guzman-Neilan bubbles."""
def __init__(self, cell, degree=None, quad_scheme=None):
def __init__(self, cell, degree=None, variant=None, quad_scheme=None):
cite("GuzmanNeilan2018")
super().__init__(FIAT.GuzmanNeilanH1div(cell, degree=degree, quad_scheme=quad_scheme))
super().__init__(FIAT.GuzmanNeilanH1div(cell, degree=degree, variant=variant, quad_scheme=quad_scheme))
Loading