diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py index c211ccc4a..c044a5922 100644 --- a/FIAT/crouzeix_raviart.py +++ b/FIAT/crouzeix_raviart.py @@ -10,109 +10,51 @@ # Last changed: 2010-01-28 from FIAT import finite_element, polynomial_set, dual_set, functional -from FIAT.check_format_variant import parse_lagrange_variant - - -def _initialize_entity_ids(topology): - entity_ids = {} - for (i, entity) in list(topology.items()): - entity_ids[i] = {} - for j in entity: - entity_ids[i][j] = [] - return entity_ids class CrouzeixRaviartDualSet(dual_set.DualSet): """Dual basis for Crouzeix-Raviart element (linears continuous at boundary midpoints).""" - def __init__(self, cell, degree): + def __init__(self, cell, degree, variant=None): # Get topology dictionary - d = cell.get_spatial_dimension() - topology = cell.get_topology() + sd = cell.get_spatial_dimension() + top = cell.get_topology() # Initialize empty nodes and entity_ids - entity_ids = _initialize_entity_ids(topology) - nodes = [None for i in list(topology[d - 1].keys())] - + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] # Construct nodes and entity_ids - for i in topology[d - 1]: - - # Construct midpoint - x = cell.make_points(d - 1, i, d)[0] - + for i in top[sd - 1]: + cur = len(nodes) + pts = cell.make_points(sd-1, i, degree+sd-1, variant=variant) # Degree of freedom number i is evaluation at midpoint - nodes[i] = functional.PointEvaluation(cell, x) - entity_ids[d - 1][i] += [i] - - # Initialize super-class - super(CrouzeixRaviartDualSet, self).__init__(nodes, cell, entity_ids) - -class CrouzeixRaviartThreeDualSet(dual_set.DualSet): - """Dual basis for Crouzeix-Raviart element (linears continuous at - boundary midpoints).""" - - def __init__(self, cell, degree): + nodes.extend(functional.PointEvaluation(cell, x) for x in pts) + entity_ids[sd - 1][i].extend(range(cur, len(nodes))) - # Get topology dictionary - d = cell.get_spatial_dimension() - topology = cell.get_topology() + cur = len(nodes) + pts = cell.make_points(sd, 0, degree, variant=variant) + nodes.extend(functional.PointEvaluation(cell, x) for x in pts) + entity_ids[sd][0].extend(range(cur, len(nodes))) - # Initialize empty nodes and entity_ids - entity_ids = _initialize_entity_ids(topology) - nodes = [None for i in range(10)] - # Construct nodes and entity_ids - entity_permutations = {0: {0: {0: []}, 1: {0: []}, 2: {0: []}}, - 1: {0: {0: [0, 1, 2], 1: [2, 1, 0]}, 1: {0: [0, 1, 2], 1: [2, 1, 0]}, 2: {0: [0, 1, 2], 1: [2, 1, 0]}}, - 2: {0: {0: [0], 1: [0], 2: [0], 3: [0], 4: [0], 5: [0]}} - } - dof_counter = 0 - for i in topology[d - 1]: - x = cell.make_points(d-1,i,4, variant="equi") - # Degree of freedom number i is evaluation at midpoint - nodes[dof_counter] = functional.PointEvaluation(cell, x[1]) - nodes[dof_counter+1] = functional.PointEvaluation(cell, x[0]) - nodes[dof_counter+2] = functional.PointEvaluation(cell, x[2]) - entity_ids[d - 1][i] += [dof_counter] - entity_ids[d - 1][i] += [dof_counter+1] - entity_ids[d - 1][i] += [dof_counter+2] - dof_counter += 3 - x = cell.make_points(d,0,3, variant="equi") - nodes[-1] = functional.PointEvaluation(cell,x[0]) - entity_ids[2][0] += [9] # Initialize super-class - super(CrouzeixRaviartThreeDualSet, self).__init__(nodes, cell, entity_ids, entity_permutations) - + super().__init__(nodes, cell, entity_ids) class CrouzeixRaviart(finite_element.CiarletElement): """The Crouzeix-Raviart finite element: K: Triangle/Tetrahedron - Polynomial space: P_1 + Polynomial space: P_1 or P_3 Dual basis: Evaluation at facet midpoints """ - def __init__(self, cell, degree,variant=None): - - # Crouzeix Raviart is only defined for polynomial degree == 1 - if not (degree in [1,3]): - raise Exception("Crouzeix-Raviart only defined for degree 1") - # Construct polynomial spaces, dual basis and initialize - # FiniteElement - if degree == 1: - space = polynomial_set.ONPolynomialSet(cell, 1) - dual = CrouzeixRaviartDualSet(cell, 1) - super(CrouzeixRaviart, self).__init__(space, dual, 1) - elif degree == 3: - ref_el = cell - splitting, point_variant = parse_lagrange_variant(variant) - if splitting is not None: - ref_el = splitting(ref_el) - poly_variant = "bubble" if ref_el.is_macrocell() else None - space = polynomial_set.ONPolynomialSet(ref_el, degree, variant=poly_variant) - formdegree = 0 # 0-form - dual = CrouzeixRaviartThreeDualSet(cell, degree) - super(CrouzeixRaviart, self).__init__(space, dual, 3, formdegree) - + def __init__(self, cell, degree, variant=None): + # Crouzeix Raviart is only defined for polynomial degree == 1 or 3 + if degree not in [1, 3]: + raise Exception("Crouzeix-Raviart only defined for degree 1 or 3") + # Construct polynomial spaces, dual basis and initialize FiniteElement + space = polynomial_set.ONPolynomialSet(cell, degree) + dual = CrouzeixRaviartDualSet(cell, degree, variant=variant) + super().__init__(space, dual, degree) diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index f81b8aa3f..e3c823d1e 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -131,6 +131,7 @@ def __init__(self, a, b): "DiscontinuousTaylor(S, 2)", "CrouzeixRaviart(I, 1)", "CrouzeixRaviart(T, 1)", + "CrouzeixRaviart(T, 3)", "CrouzeixRaviart(S, 1)", "RaviartThomas(I, 1)", "RaviartThomas(I, 2)",