From f96b28c23b87bd1afaff2696eab9e8eb487e1349 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 20 Nov 2024 17:06:45 +0000 Subject: [PATCH] mass-decoupling variant --- FIAT/brezzi_douglas_marini.py | 4 ++-- FIAT/check_format_variant.py | 8 +++----- FIAT/demkowicz.py | 26 +++++++++++++++++--------- FIAT/hierarchical.py | 4 ++-- FIAT/nedelec.py | 4 ++-- FIAT/nedelec_second_kind.py | 4 ++-- FIAT/raviart_thomas.py | 4 ++-- test/unit/test_hct.py | 22 ++++++++++++++++++++++ 8 files changed, 52 insertions(+), 24 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index d4181de18..5f5451a85 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -94,8 +94,8 @@ def __init__(self, ref_el, degree, variant=None): sd = ref_el.get_spatial_dimension() poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, )) - if variant == "demkowicz": - dual = demkowicz.DemkowiczDual(ref_el, degree, "HDiv") + if variant and variant.startswith("demkowicz"): + dual = demkowicz.DemkowiczDual(ref_el, degree, "HDiv", variant=variant) elif variant == "fdm": dual = demkowicz.FDMDual(ref_el, degree, "HDiv", type(self)) else: diff --git a/FIAT/check_format_variant.py b/FIAT/check_format_variant.py index 30ff58f21..e390016ef 100644 --- a/FIAT/check_format_variant.py +++ b/FIAT/check_format_variant.py @@ -30,10 +30,9 @@ def check_format_variant(variant, degree): if variant is None: variant = "integral" - match = re.match(r"^integral(?:\((\d+)\))?$", variant) + match = re.match(r"^(?:(\w+))?(?:\((\d+)\))?$", variant) if match: - variant = "integral" - extra_degree, = match.groups() + variant, extra_degree = match.groups() extra_degree = int(extra_degree) if extra_degree is not None else 0 interpolant_degree = degree + extra_degree if interpolant_degree < degree: @@ -61,7 +60,7 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False): default = "integral" if integral else "spectral" if integral: - supported_point_variants = {"integral": None, "fdm": "fdm", "demkowicz": "demkowicz"} + supported_point_variants = {"integral": None, "fdm": "fdm", "demkowicz": "demkowicz", "demkowicz-mass": "demkowicz-mass"} elif discontinuous: supported_point_variants = supported_dg_variants else: @@ -71,7 +70,6 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False): splitting = None splitting_args = tuple() point_variant = supported_point_variants[default] - for pre_opt in options: opt = pre_opt.lower() if opt in supported_splits: diff --git a/FIAT/demkowicz.py b/FIAT/demkowicz.py index 54516d5df..51782a59d 100644 --- a/FIAT/demkowicz.py +++ b/FIAT/demkowicz.py @@ -65,7 +65,7 @@ def map_duals(ref_el, dim, entity, mapping, Q_ref, Phis): class DemkowiczDual(DualSet): - def __init__(self, ref_el, degree, sobolev_space, kind=None): + def __init__(self, ref_el, degree, sobolev_space, kind=None, variant=None): nodes = [] entity_ids = {} reduced_dofs = {} @@ -76,6 +76,8 @@ def __init__(self, ref_el, degree, sobolev_space, kind=None): dual_mapping = {"HCurl": "contravariant", "HDiv": "covariant"}.get(sobolev_space, None) if kind is None: kind = 1 if formdegree == 0 else 2 + if variant is None: + variant = "demkowicz" for dim in sorted(top): entity_ids[dim] = {} @@ -91,7 +93,7 @@ def __init__(self, ref_el, degree, sobolev_space, kind=None): entity_ids[dim][entity] = list(range(cur, len(nodes))) reduced_dofs[dim] = len(nodes) else: - Q_ref, Phis, rdofs = self._reference_duals(dim, degree, formdegree, sobolev_space, kind) + Q_ref, Phis, rdofs = self._reference_duals(dim, degree, formdegree, sobolev_space, kind, variant) reduced_dofs[dim] = rdofs mapping = dual_mapping if dim == sd else trace for entity in sorted(top[dim]): @@ -101,9 +103,9 @@ def __init__(self, ref_el, degree, sobolev_space, kind=None): entity_ids[dim][entity] = list(range(cur, len(nodes))) self._reduced_dofs = reduced_dofs - super(DemkowiczDual, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) - def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind): + def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind, variant): facet = symmetric_simplex(dim) Q = create_quadrature(facet, 2 * degree) Qpts, Qwts = Q.get_points(), Q.get_weights() @@ -122,10 +124,15 @@ def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind): if formdegree >= dim: K = inner(trial[:1], trial, Qwts) else: - dtrial = exterior_derivative(P_at_qpts) - if dim == 2 and formdegree == 1 and sobolev_space == "HDiv": - dtrial = dtrial[:, None, :] - K = self._bubble_derivative_moments(facet, degree, formdegree, kind, Qpts, Qwts, dtrial) + if variant == "demkowicz": + deriv = True + dtrial = exterior_derivative(P_at_qpts) + if dim == 2 and formdegree == 1 and sobolev_space == "HDiv": + dtrial = dtrial[:, None, :] + else: + deriv = False + dtrial = trial + K = self._bubble_derivative_moments(facet, degree, formdegree, kind, Qpts, Qwts, dtrial, deriv=deriv) reduced_dofs = K.shape[0] # Evaluate type-II degrees of freedom on P @@ -220,7 +227,7 @@ def __init__(self, ref_el, degree, sobolev_space, element): self.Q = Q self.V0 = phis[(0,) * sd] self.V1 = exterior_derivative(phis) - super(FDMDual, self).__init__(ref_el, degree, sobolev_space, kind=None) + super().__init__(ref_el, degree, sobolev_space, kind=None) def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind): entity_dofs = self.fe.entity_dofs() @@ -286,6 +293,7 @@ def _reference_duals(self, dim, degree, formdegree, sobolev_space, kind): kind = 1 variant = "fdm" variant = "demkowicz" + variant = "demkowicz-mass" # variant = None space_dict = {"H1": (CG, grad), "HCurl": (N1Curl if kind == 1 else N2Curl, curl), diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py index 945acde6c..f129839e3 100644 --- a/FIAT/hierarchical.py +++ b/FIAT/hierarchical.py @@ -78,7 +78,7 @@ def __init__(self, ref_el, degree, variant=None): if splitting is not None: ref_el = splitting(ref_el) poly_set = ONPolynomialSet(ref_el, degree) - if variant == "demkowicz": + if variant and variant.startswith("demkowicz"): dual = demkowicz.DemkowiczDual(ref_el, degree, "L2") elif variant == "fdm": dual = demkowicz.FDMDual(ref_el, degree, "L2", type(self)) @@ -135,7 +135,7 @@ def __init__(self, ref_el, degree, variant=None): if degree < 1: raise ValueError(f"{type(self).__name__} elements only valid for k >= 1") poly_set = ONPolynomialSet(ref_el, degree, variant="bubble") - if variant == "demkowicz": + if variant and variant.startswith("demkowicz"): dual = demkowicz.DemkowiczDual(ref_el, degree, "H1") elif variant == "fdm": dual = demkowicz.FDMDual(ref_el, degree, "H1", type(self)) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index d61712c6d..6471eca58 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -196,8 +196,8 @@ class Nedelec(finite_element.CiarletElement): def __init__(self, ref_el, degree, variant=None): - if variant == "demkowicz": - dual = demkowicz.DemkowiczDual(ref_el, degree, "HCurl", kind=1) + if variant and variant.startswith("demkowicz"): + dual = demkowicz.DemkowiczDual(ref_el, degree, "HCurl", kind=1, variant=variant) elif variant == "fdm": dual = demkowicz.FDMDual(ref_el, degree, "HCurl", type(self)) else: diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index c43d05969..0b6e92618 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -198,8 +198,8 @@ def __init__(self, ref_el, degree, variant=None): sd = ref_el.get_spatial_dimension() poly_set = ONPolynomialSet(ref_el, degree, (sd, ), variant="bubble") - if variant == "demkowicz": - dual = demkowicz.DemkowiczDual(ref_el, degree, "HCurl") + if variant and variant.startswith("demkowicz"): + dual = demkowicz.DemkowiczDual(ref_el, degree, "HCurl", variant=variant) elif variant == "fdm": dual = demkowicz.FDMDual(ref_el, degree, "HCurl", type(self)) else: diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index a450e3326..2438cdca9 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -141,8 +141,8 @@ class RaviartThomas(finite_element.CiarletElement): def __init__(self, ref_el, degree, variant=None): - if variant == "demkowicz": - dual = demkowicz.DemkowiczDual(ref_el, degree, "HDiv", kind=1) + if variant and variant.startswith("demkowicz"): + dual = demkowicz.DemkowiczDual(ref_el, degree, "HDiv", kind=1, variant=variant) elif variant == "fdm": dual = demkowicz.FDMDual(ref_el, degree, "HDiv", type(self)) else: diff --git a/test/unit/test_hct.py b/test/unit/test_hct.py index 84e0a2379..7e462a65c 100644 --- a/test/unit/test_hct.py +++ b/test/unit/test_hct.py @@ -5,6 +5,8 @@ from FIAT.reference_element import ufc_simplex from FIAT.functional import PointEvaluation from FIAT.macro import CkPolynomialSet +from FIAT.quadrature_schemes import create_quadrature +from FIAT.jacobi import eval_jacobi @pytest.fixture @@ -74,3 +76,23 @@ def test_full_polynomials(cell, reduced): C1 = CkPolynomialSet(ref_complex, degree, order=1, variant="bubble") C1_tab = C1.tabulate(pts)[(0, 0)] assert span_greater_equal(tab, C1_tab) + + +def test_reduced_normal_derivative(cell): + fe = HCT(cell, reduced=True) + + ref_line = cell.construct_subelement(1) + Q = create_quadrature(ref_line, fe.degree()+2) + qpts, qwts = Q.get_points(), Q.get_weights() + + bary = ref_line.compute_barycentric_coordinates(qpts) + leg2 = eval_jacobi(0, 0, 2, bary[:, 1] - bary[:, 0]) + wts = numpy.multiply(leg2, qwts) + top = cell.get_topology() + for e in top[1]: + n = cell.compute_normal(e) + vals = fe.tabulate(1, qpts, entity=(1, e)) + fn = vals[(1, 0)] * n[0] + vals[(0, 1)] * n[1] + + fn = fn[:-3] + assert numpy.allclose(numpy.dot(fn, wts), 0)