Skip to content

Symmetric tangential constraints #248

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

Open
pbrubeck wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Symmetric tangential constraints #248

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
19 changes: 15 additions & 4 deletions FIAT/bernardi_raugel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,15 +74,26 @@ def __init__(self, ref_el, order=1, degree=None, reduced=False, ref_complex=None
# Face moments of tangential components against dual bubbles
ref_facet = ref_complex.construct_subcomplex(sd-1)
ref_area = ref_facet.volume()

# Quadrature and weight function for tangential constraints
codim = sd-1 if degree == 1 and ref_facet.is_macrocell() else 0
Qt_ref, phis = make_dual_bubbles(ref_facet, degree, codim=codim, scale=1)
ft_at_qpts = phis[-1]
if codim == 0:
scale = (-1)**(sd-1) * 0.5 * ref_area
Qt_ref = parse_quadrature_scheme(ref_facet, 2*degree, quad_scheme=None)
Pk = polynomial_set.ONPolynomialSet(ref_facet, degree, scale="orthonormal")
coeffs = polynomial_set.project(ref_facet.compute_bubble, Pk, Qt_ref)
coeffs[..., :expansions.polynomial_dimension(ref_facet, degree-1)] = 0
Pt = polynomial_set.PolynomialSet(ref_facet, degree, degree,
Pk.get_expansion_set(),
coeffs)
ft_at_qpts = Pt.tabulate(Qt_ref.get_points())[(0,)*(sd-1)][-1]
bf_at_qpts = ref_facet.compute_bubble(Qt_ref.get_points())
bf_wts = numpy.multiply(bf_at_qpts, Qt_ref.get_weights())
scale = numpy.sum(bf_wts) / numpy.dot(ft_at_qpts, bf_wts)
scale /= ref_area
else:
Qt_ref, phis = make_dual_bubbles(ref_facet, degree, codim=codim, scale=1)
ft_at_qpts = phis[-1]
scale = ref_area / numpy.dot(ft_at_qpts, Qt_ref.get_weights())

ft_at_qpts *= scale
if codim != 0:
ft_at_qpts -= numpy.dot(ft_at_qpts, Qt_ref.get_weights()) / ref_area
Expand Down
9 changes: 6 additions & 3 deletions FIAT/polynomial_set.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,11 +139,14 @@ def project(f, U, Q):
a polynomial set U. Numerical integration is performed by
quadrature rule Q.
"""
sd = Q.ref_el.get_spatial_dimension()
pts = Q.get_points()
wts = Q.get_weights()
f_at_qps = [f(x) for x in pts]
U_at_qps = U.tabulate(pts)
coeffs = numpy.array([sum(wts * f_at_qps * phi) for phi in U_at_qps])
U_at_qps = U.tabulate(pts)[(0,)*sd]
f_at_qps = numpy.reshape(f(pts), (-1, *U_at_qps.shape[1:]))
duals = numpy.multiply(U_at_qps, wts)
coeffs = numpy.tensordot(f_at_qps, duals,
axes=(range(1, duals.ndim),)*2)
return coeffs


Expand Down
Loading