crouzeix_raviart.py

# Copyright (C) 2010 Marie E. Rognes
#
# This file is part of FIAT (https://www.fenicsproject.org)
#
# SPDX-License-Identifier:    LGPL-3.0-or-later
#
# Written by Marie E. Rognes <meg@simula.no> based on original
# implementation by Robert C. Kirby.
#
# Last changed: 2010-01-28

import numpy
from FIAT import finite_element, polynomial_set, dual_set, functional, macro
from FIAT.check_format_variant import check_format_variant, parse_quadrature_scheme
from FIAT.quadrature import FacetQuadratureRule


class CrouzeixRaviartDualSet(dual_set.DualSet):

    def __init__(self, ref_el, degree, variant, interpolant_deg, quad_scheme):
        # Get topology dictionary
        sd = ref_el.get_spatial_dimension()
        top = ref_el.get_topology()

        if degree > 1 and sd != 2:
            raise NotImplementedError("High-order Crouzeix-Raviart is only implemented on triangles.")

        # Initialize empty nodes and entity_ids
        entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top}
        nodes = []

        # Construct nodes and entity_ids
        if variant == "integral":
            for dim in sorted(top):
                if dim == 0 and dim != sd-1:
                    # Skip vertex dofs
                    continue
                facet = ref_el.construct_subelement(dim)
                if dim == 0:
                    Q_facet = parse_quadrature_scheme(facet, degree + interpolant_deg-1, quad_scheme)
                    phis = numpy.ones((1, len(Q_facet.pts)))
                else:
                    k = degree - 1 if dim == sd-1 else degree - (1+dim)
                    if k < 0:
                        continue
                    Q_facet = parse_quadrature_scheme(facet, k + interpolant_deg, quad_scheme)
                    poly_set = polynomial_set.ONPolynomialSet(facet, k)
                    phis = poly_set.tabulate(Q_facet.get_points())[(0,) * dim]

                for i in sorted(top[dim]):
                    cur = len(nodes)
                    Q = FacetQuadratureRule(ref_el, dim, i, Q_facet, avg=True)
                    nodes.extend(functional.IntegralMoment(ref_el, Q, phi) for phi in phis)
                    entity_ids[dim][i].extend(range(cur, len(nodes)))
        else:
            for dim in sorted(top):
                if dim == 0 and dim != sd-1:
                    # Skip vertex dofs
                    continue
                for i in sorted(top[dim]):
                    cur = len(nodes)
                    if dim == sd-1 and dim != 0:
                        pts = ref_el.make_points(dim, i, degree-1, variant="gl", interior=0)
                    else:
                        pts = ref_el.make_points(dim, i, degree, variant="gll")
                    nodes.extend(functional.PointEvaluation(ref_el, x) for x in pts)
                    entity_ids[dim][i].extend(range(cur, len(nodes)))

        super().__init__(nodes, ref_el, entity_ids)


class CrouzeixRaviart(finite_element.CiarletElement):
    """The Crouzeix-Raviart finite element:

    K:                 Triangle/Tetrahedron
    Polynomial space:  P_k
    Dual basis:        Evaluation at points or integral moments
    """
    DEFAULT_DEGREE = 1

    def __init__(self, ref_el, degree, variant=None, quad_scheme=None):
        degree = self._parse_degree(degree)
        if degree % 2 != 1:
            raise ValueError("Crouzeix-Raviart only defined for odd degree")

        splitting, variant, interpolant_deg = check_format_variant(variant, degree)
        if splitting is not None:
            ref_el = splitting(ref_el)

        if ref_el.is_macrocell():
            base_element = type(self)(ref_el.get_parent(), degree)
            poly_set = macro.MacroPolynomialSet(ref_el, base_element)
        else:
            poly_set = polynomial_set.ONPolynomialSet(ref_el, degree)

        dual = CrouzeixRaviartDualSet(ref_el, degree, variant, interpolant_deg, quad_scheme)
        super().__init__(poly_set, dual, degree)