diff --git a/FIAT/brezzi_douglas_fortin_marini.py b/FIAT/brezzi_douglas_fortin_marini.py index 4361e70f6..b6f8526ec 100644 --- a/FIAT/brezzi_douglas_fortin_marini.py +++ b/FIAT/brezzi_douglas_fortin_marini.py @@ -1,116 +1,25 @@ -from FIAT import (finite_element, functional, dual_set, - polynomial_set, lagrange) +from FIAT.expansions import polynomial_dimension +from FIAT.brezzi_douglas_marini import BrezziDouglasMarini +from FIAT.nodal_enriched import NodalEnrichedElement +from FIAT.restricted import RestrictedElement -import numpy - -class BDFMDualSet(dual_set.DualSet): - def __init__(self, ref_el, degree): - - # Initialize containers for map: mesh_entity -> dof number and - # dual basis - entity_ids = {} - nodes = [] - - sd = ref_el.get_spatial_dimension() - t = ref_el.get_topology() - - # Define each functional for the dual set - # codimension 1 facet normals. - # note this will die for degree greater than 1. - for i in range(len(t[sd - 1])): - pts_cur = ref_el.make_points(sd - 1, i, sd + degree) - for j in range(len(pts_cur)): - pt_cur = pts_cur[j] - f = functional.PointScaledNormalEvaluation(ref_el, i, pt_cur) - nodes.append(f) - - # codimension 1 facet tangents. - # because the tangent component is discontinuous, these actually - # count as internal nodes. - tangent_count = 0 - for i in range(len(t[sd - 1])): - pts_cur = ref_el.make_points(sd - 1, i, sd + degree - 1) - tangent_count += len(pts_cur) - for j in range(len(pts_cur)): - pt_cur = pts_cur[j] - f = functional.PointEdgeTangentEvaluation(ref_el, i, pt_cur) - nodes.append(f) - - # sets vertices (and in 3d, edges) to have no nodes - for i in range(sd - 1): - entity_ids[i] = {} - for j in range(len(t[i])): - entity_ids[i][j] = [] - - cur = 0 - - # set codimension 1 (edges 2d, faces 3d) dof - pts_facet_0 = ref_el.make_points(sd - 1, 0, sd + degree) - pts_per_facet = len(pts_facet_0) - - entity_ids[sd - 1] = {} - for i in range(len(t[sd - 1])): - entity_ids[sd - 1][i] = list(range(cur, cur + pts_per_facet)) - cur += pts_per_facet - - # internal nodes - entity_ids[sd] = {0: list(range(cur, cur + tangent_count))} - cur += tangent_count - - super().__init__(nodes, ref_el, entity_ids) - - -def BDFMSpace(ref_el, order): - sd = ref_el.get_spatial_dimension() - if sd != 2: - raise Exception("BDFM_k elements only valid for dim 2") - # Note that order will be 2. - - # Linear vector valued space. Since the embedding degree of this element - # is 2, this is implemented by taking the quadratic space and selecting - # the linear polynomials. - vec_poly_set = polynomial_set.ONPolynomialSet(ref_el, order, (sd,)) - # Linears are the first three polynomials in each dimension. - vec_poly_set = vec_poly_set.take([0, 1, 2, 6, 7, 8]) - - # Scalar quadratic Lagrange element. - lagrange_ele = lagrange.Lagrange(ref_el, order) - # Select the dofs associated with the edges. - edge_dofs_dict = lagrange_ele.dual.get_entity_ids()[sd - 1] - edge_dofs = numpy.array([(edge, dof) - for edge, dofs in list(edge_dofs_dict.items()) - for dof in dofs]) - - tangent_polys = lagrange_ele.poly_set.take(edge_dofs[:, 1]) - new_coeffs = numpy.zeros((tangent_polys.get_num_members(), sd, tangent_polys.coeffs.shape[-1])) - - # Outer product of the tangent vectors with the quadratic edge polynomials. - for i, (edge, dof) in enumerate(edge_dofs): - tangent = ref_el.compute_edge_tangent(edge) - - new_coeffs[i, :, :] = numpy.outer(tangent, tangent_polys.coeffs[i, :]) - - bubble_set = polynomial_set.PolynomialSet(ref_el, - order, - order, - vec_poly_set.get_expansion_set(), - new_coeffs) - - element_set = polynomial_set.polynomial_set_union_normalized(bubble_set, vec_poly_set) - return element_set - - -class BrezziDouglasFortinMarini(finite_element.CiarletElement): +def BrezziDouglasFortinMarini(ref_el, degree, variant=None, quad_scheme=None): """The BDFM element""" - - def __init__(self, ref_el, degree): - - if degree != 2: - raise Exception("BDFM_k elements only valid for k == 2") - - poly_set = BDFMSpace(ref_el, degree) - dual = BDFMDualSet(ref_el, degree - 1) - formdegree = ref_el.get_spatial_dimension() - 1 - super().__init__(poly_set, dual, degree, formdegree, - mapping="contravariant piola") + if variant == "point": + BDM_I = RestrictedElement(BrezziDouglasMarini(ref_el, degree, variant=variant), restriction_domain="interior") + BDM_F = RestrictedElement(BrezziDouglasMarini(ref_el, degree-1, variant=variant), restriction_domain="facet") + return NodalEnrichedElement(BDM_I, BDM_F) + else: + BDM = BrezziDouglasMarini(ref_el, degree, variant=variant, quad_scheme=quad_scheme) + entity_ids = BDM.dual.get_entity_ids() + sd = ref_el.get_spatial_dimension() + indices = [] + for dim in sorted(entity_ids): + if dim == sd-1: + s = slice(polynomial_dimension(ref_el.construct_subelement(dim), degree-1)) + else: + s = slice(None) + for entity in sorted(entity_ids[dim]): + indices.extend(entity_ids[dim][entity][s]) + return RestrictedElement(BDM, indices) diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py index 3ce7fe2d2..6498298b1 100644 --- a/finat/fiat_elements.py +++ b/finat/fiat_elements.py @@ -460,8 +460,8 @@ def entity_permutations(self): class BrezziDouglasFortinMarini(VectorFiatElement): - def __init__(self, cell, degree): - super().__init__(FIAT.BrezziDouglasFortinMarini(cell, degree)) + def __init__(self, cell, degree, **kwargs): + super().__init__(FIAT.BrezziDouglasFortinMarini(cell, degree, **kwargs)) class Nedelec(VectorFiatElement): diff --git a/test/FIAT/regression/fiat-reference-data-id b/test/FIAT/regression/fiat-reference-data-id index 6236063fa..dfbfe6615 100644 --- a/test/FIAT/regression/fiat-reference-data-id +++ b/test/FIAT/regression/fiat-reference-data-id @@ -1 +1 @@ -0c8c97f7e4919402129e5ff3b54e3f0b9e902b7c +508bd755e024010f6fc691a36e51a8f4d7de7efe diff --git a/test/FIAT/regression/scripts/getdata b/test/FIAT/regression/scripts/getdata index 9eacdac3f..cebdb5945 100755 --- a/test/FIAT/regression/scripts/getdata +++ b/test/FIAT/regression/scripts/getdata @@ -30,5 +30,5 @@ source scripts/parameters DATA_ID=$1 && [ -z "$DATA_ID" ] && DATA_ID=`cat $DATA_ID_FILE` # Checkout data referenced by id -(cd $DATA_DIR && git checkout -B auto $DATA_ID) +(cd $DATA_DIR && git checkout -B main $DATA_ID) exit $? diff --git a/test/FIAT/regression/scripts/parameters b/test/FIAT/regression/scripts/parameters index b8b0ce849..857721506 100755 --- a/test/FIAT/regression/scripts/parameters +++ b/test/FIAT/regression/scripts/parameters @@ -1,6 +1,6 @@ #OUTPUT_DIR="output" -[ ! -z ${DATA_REPO_GIT+x} ] || DATA_REPO_GIT="git@bitbucket.org:fenics-project/fiat-reference-data.git" -DATA_REPO_HTTPS="https://bitbucket.org/fenics-project/fiat-reference-data.git" +[ ! -z ${DATA_REPO_GIT+x} ] || DATA_REPO_GIT="git@github.com:firedrakeproject/fiat-reference-data.git" +DATA_REPO_HTTPS="https://github.com/firedrakeproject/fiat-reference-data.git" DATA_DIR="fiat-reference-data" DATA_ID_FILE="fiat-reference-data-id" OUTPUT_DIR="$DATA_DIR" diff --git a/test/FIAT/regression/test_regression.py b/test/FIAT/regression/test_regression.py index fc5ee01e5..4af77dfd8 100644 --- a/test/FIAT/regression/test_regression.py +++ b/test/FIAT/regression/test_regression.py @@ -224,7 +224,10 @@ def quadrature_test_case_name(test_case): ("Brezzi-Douglas-Marini", 3, 1), ("Brezzi-Douglas-Marini", 3, 2), ("Brezzi-Douglas-Marini", 3, 3), + ("Brezzi-Douglas-Fortin-Marini", 2, 1), ("Brezzi-Douglas-Fortin-Marini", 2, 2), + ("Brezzi-Douglas-Fortin-Marini", 3, 1), + ("Brezzi-Douglas-Fortin-Marini", 3, 2), ("Raviart-Thomas", 2, 1), ("Raviart-Thomas", 2, 2), ("Raviart-Thomas", 2, 3), diff --git a/test/FIAT/unit/test_bdfm.py b/test/FIAT/unit/test_bdfm.py new file mode 100644 index 000000000..ef90fdc36 --- /dev/null +++ b/test/FIAT/unit/test_bdfm.py @@ -0,0 +1,55 @@ +import pytest +import numpy +from FIAT import BrezziDouglasFortinMarini +from FIAT.reference_element import ufc_simplex +from FIAT.polynomial_set import ONPolynomialSet +from FIAT.expansions import polynomial_dimension +from FIAT.quadrature_schemes import create_quadrature +from FIAT.quadrature import FacetQuadratureRule + + +@pytest.fixture(params=[2, 3]) +def cell(request): + dim = request.param + K = ufc_simplex(dim) + return K + + +def inner(u, v, qwts): + return numpy.tensordot(numpy.multiply(u, qwts), v, axes=(tuple(range(1, u.ndim)), )*2) + + +@pytest.mark.parametrize("degree", (2, 4)) +@pytest.mark.parametrize("variant", ("point", "integral", "integral(1)")) +def test_bdfm_space(cell, degree, variant): + + fe = BrezziDouglasFortinMarini(cell, degree, variant=variant) + + sd = cell.get_spatial_dimension() + ref_face = cell.construct_subelement(sd-1) + Q_face = create_quadrature(ref_face, 2*degree) + + dimPk = polynomial_dimension(ref_face, degree) + dimPkm1 = polynomial_dimension(ref_face, degree-1) + Phis = ONPolynomialSet(ref_face, degree) + Phis = Phis.take(range(dimPkm1, dimPk)) + phis_at_qpts = Phis.tabulate(Q_face.get_points())[(0,)*(sd-1)] + + expected_dim = sd*polynomial_dimension(cell, degree) - (sd+1)*len(Phis) + assert fe.space_dimension() == expected_dim + + for f in cell.topology[sd-1]: + Q = FacetQuadratureRule(cell, sd-1, f, Q_face) + vals = fe.tabulate(0, Q.get_points())[(0,)*sd] + + # Assert that the normal component is one degree lower + n = cell.compute_normal(f) + normal_trace = numpy.tensordot(vals, n, axes=(1, 0)) + normal_moments = inner(normal_trace, phis_at_qpts, Q.get_weights()) + assert numpy.allclose(normal_moments, 0) + + # Assert that the tangential component has the expected degree + for t in cell.compute_tangents(sd-1, f): + tangential_trace = numpy.tensordot(vals, t, axes=(1, 0)) + tangential_moments = inner(tangential_trace, phis_at_qpts, Q.get_weights()) + assert not numpy.allclose(tangential_moments, 0) diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py index 7e2586097..361cd510e 100644 --- a/test/FIAT/unit/test_fiat.py +++ b/test/FIAT/unit/test_fiat.py @@ -300,6 +300,9 @@ def __init__(self, a, b): "GopalakrishnanLedererSchoberlSecondKind(S, 1)", "GopalakrishnanLedererSchoberlSecondKind(S, 2)", "BrezziDouglasFortinMarini(T, 2)", + "BrezziDouglasFortinMarini(T, 3)", + "BrezziDouglasFortinMarini(S, 2)", + "BrezziDouglasFortinMarini(T, 2, variant='point')", "GaussLegendre(I, 0)", "GaussLegendre(I, 1)", "GaussLegendre(I, 2)", diff --git a/test/finat/test_create_fiat_element.py b/test/finat/test_create_fiat_element.py index 00b53b24f..64490d957 100644 --- a/test/finat/test_create_fiat_element.py +++ b/test/finat/test_create_fiat_element.py @@ -11,7 +11,6 @@ supported_elements = { # These all map directly to FIAT elements "Brezzi-Douglas-Marini": FIAT.BrezziDouglasMarini, - "Brezzi-Douglas-Fortin-Marini": FIAT.BrezziDouglasFortinMarini, "Lagrange": FIAT.Lagrange, "Nedelec 1st kind H(curl)": FIAT.Nedelec, "Nedelec 2nd kind H(curl)": FIAT.NedelecSecondKind, @@ -29,7 +28,6 @@ def create_element(ufl_element): @pytest.fixture(params=["BDM", - "BDFM", "Lagrange", "N1curl", "N2curl",