diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 2c9aeb89..a3f1fbb6 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -30,13 +30,12 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # Facet nodes are \int_F v\cdot n p ds where p \in P_{q} # degree is q Q_ref = create_quadrature(facet, interpolant_deg + degree) - Pq = polynomial_set.ONPolynomialSet(facet, degree) + Pq = polynomial_set.ONPolynomialSet(facet, degree, scale=1) Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] for f in top[sd - 1]: cur = len(nodes) Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref) - Jdet = Q.jacobian_determinant() - n = ref_el.compute_scaled_normal(f) / Jdet + n = ref_el.compute_normal(f) phis = n[None, :, None] * Pq_at_qpts[:, None, :] nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 854855d8..acaf14e8 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -252,8 +252,8 @@ def __init__(self, ref_el, scale=None, variant=None): self.mapping = lambda x: numpy.dot(self.A, x) + self.b self._dmats_cache = {} if scale is None: - scale = math.sqrt(1.0 / self.base_ref_el.volume()) - elif isinstance(scale, str): + scale = "orthonormal" + if isinstance(scale, str): scale = scale.lower() if scale == "orthonormal": scale = math.sqrt(1.0 / ref_el.volume()) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index c0513302..c36e2b15 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -123,7 +123,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): if phi_deg >= 0: facet = ref_el.construct_subelement(dim) Q_ref = create_quadrature(facet, interpolant_deg + phi_deg) - Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,)) + Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,), scale="L2 piola") Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * dim] Phis = numpy.transpose(Phis, (0, 2, 1)) @@ -165,7 +165,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): interpolant_deg = degree cur = len(nodes) Q = create_quadrature(ref_el, interpolant_deg + phi_deg) - Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg) + Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg, scale="L2 piola") Phis = Pqmd.tabulate(Q.get_points())[(0,) * dim] nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (dim,)) for d in range(dim) for phi in Phis) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 6e526972..08e5c137 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -129,7 +129,7 @@ def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant ref_facet = cell.construct_subelement(codim) Q_ref = create_quadrature(ref_facet, interpolant_deg + rt_degree) if codim == 1: - Phi = ONPolynomialSet(ref_facet, rt_degree, (codim,)) + Phi = ONPolynomialSet(ref_facet, rt_degree, (codim,), scale="L2 piola") else: # Construct Raviart-Thomas on the reference facet RT = RaviartThomas(ref_facet, rt_degree, variant) @@ -149,10 +149,10 @@ def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant # Get the quadrature and Jacobian on this facet Q_facet = FacetQuadratureRule(cell, codim, facet, Q_ref) J = Q_facet.jacobian() - detJ = Q_facet.jacobian_determinant() + Jdet = Q_facet.jacobian_determinant() # Map Phis -> phis (reference values to physical values) - piola_map = J / detJ + piola_map = J / Jdet phis = numpy.dot(Phis, piola_map.T) phis = numpy.transpose(phis, (0, 2, 1)) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 64271edd..dd97c486 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -122,7 +122,6 @@ class ONPolynomialSet(PolynomialSet): """ def __init__(self, ref_el, degree, shape=tuple(), scale=None, variant=None): - if shape == tuple(): num_components = 1 else: diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index c66c4ff4..1f190ce5 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -75,13 +75,12 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1} # degree is q - 1 Q_ref = create_quadrature(facet, interpolant_deg + degree - 1) - Pq = polynomial_set.ONPolynomialSet(facet, degree - 1) + Pq = polynomial_set.ONPolynomialSet(facet, degree - 1, scale=1) Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] for f in top[sd - 1]: cur = len(nodes) Q = FacetQuadratureRule(ref_el, sd-1, f, Q_ref) - Jdet = Q.jacobian_determinant() - n = ref_el.compute_scaled_normal(f) / Jdet + n = ref_el.compute_normal(f) phis = n[None, :, None] * Pq_at_qpts[:, None, :] nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) @@ -91,7 +90,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): if degree > 1: cur = len(nodes) Q = create_quadrature(ref_el, interpolant_deg + degree - 2) - Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 2) + Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 2, scale=1) Pkm1_at_qpts = Pkm1.tabulate(Q.get_points())[(0,) * sd] nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) for d in range(sd) diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 5f0b0f28..13d5652e 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -530,20 +530,22 @@ def test_error_point_high_order(element): eval(element) -@pytest.mark.parametrize('cell', [I, T, S]) -def test_expansion_orthonormality(cell): - from FIAT import expansions +@pytest.mark.parametrize('degree', [0, 10]) +@pytest.mark.parametrize('dim', range(1, 4)) +@pytest.mark.parametrize('cell', ["default", "ufc"]) +def test_expansion_orthonormality(cell, dim, degree): + from FIAT.expansions import ExpansionSet from FIAT.quadrature_schemes import create_quadrature - U = expansions.ExpansionSet(cell) - degree = 10 - rule = create_quadrature(cell, 2*degree) - phi = U.tabulate(degree, rule.pts) - qwts = rule.get_weights() - scale = 2 ** cell.get_spatial_dimension() - results = scale * np.dot(np.multiply(phi, qwts), phi.T) - - assert np.allclose(results, np.diag(np.diag(results))) - assert np.allclose(np.diag(results), 1.0) + from FIAT import reference_element + make_cell = {"default": reference_element.default_simplex, + "ufc": reference_element.ufc_simplex}[cell] + ref_el = make_cell(dim) + U = ExpansionSet(ref_el) + Q = create_quadrature(ref_el, 2 * degree) + qpts, qwts = Q.get_points(), Q.get_weights() + phi = U.tabulate(degree, qpts) + results = np.dot(np.multiply(phi, qwts), phi.T) + assert np.allclose(results, np.eye(results.shape[0])) @pytest.mark.parametrize('dim', range(1, 4))