diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 04ac6f3a..8fb54397 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -14,7 +14,7 @@ from FIAT.alfeld_sorokina import AlfeldSorokina from FIAT.arnold_qin import ArnoldQin from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanH1div -from FIAT.stokes import Stokes, MacroStokes +from FIAT.stokes import Stokes, MacroStokes, DivStokes from FIAT.christiansen_hu import ChristiansenHu from FIAT.johnson_mercier import JohnsonMercier from FIAT.brezzi_douglas_marini import BrezziDouglasMarini @@ -107,6 +107,7 @@ "Regge": Regge, "Stokes": Stokes, "Macro Stokes": MacroStokes, + "Div Stokes": DivStokes, "EnrichedElement": EnrichedElement, "NodalEnrichedElement": NodalEnrichedElement, "QuadraticPowellSabin6": QuadraticPowellSabin6, diff --git a/FIAT/stokes.py b/FIAT/stokes.py index d5e40f7e..7a0bf9f3 100644 --- a/FIAT/stokes.py +++ b/FIAT/stokes.py @@ -8,7 +8,8 @@ import scipy from FIAT import finite_element, dual_set, expansions, jacobi, macro -from FIAT.functional import (ComponentPointEvaluation, +from FIAT.functional import (PointEvaluation, + ComponentPointEvaluation, PointTangentialDerivative, PointTangentialSecondDerivative, PointSecondDerivative, @@ -18,6 +19,7 @@ from FIAT.polynomial_set import make_bubbles, ONPolynomialSet from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature +from FIAT.bernstein import Bernstein def eps(table_u): @@ -42,13 +44,24 @@ def map_duals(ref_el, dim, entity, mapping, Q_ref, Phis): return Q, phis -def jacobi_duals(ref_el, weight, moment_degree, degree): +def jacobi_duals(ref_el, weight, trial_degree, test_degree): facet = ref_el.construct_subelement(1) - Q = create_quadrature(facet, degree + moment_degree) + Q = create_quadrature(facet, trial_degree + test_degree) x = facet.compute_barycentric_coordinates(Q.get_points()) xhat = x[:, 1:] - x[:, :1] a = weight - phis = jacobi.eval_jacobi_batch(a, a, moment_degree, xhat) + phis = jacobi.eval_jacobi_batch(a, a, test_degree, xhat) + return Q, phis + + +def dubiner_duals(ref_el, dim, trial_degree, test_degree): + if dim == 0: + return None, [] + + facet = ref_el.construct_subelement(dim) + V = ONPolynomialSet(facet, test_degree, scale="orthonormal") + Q = create_quadrature(facet, trial_degree + V.degree) + phis = V.tabulate(Q.get_points())[(0,)*dim] return Q, phis @@ -64,7 +77,7 @@ def __init__(self, ref_el, degree): moment_degree = degree - 2*sd edge_weight = sd-1 - Q_edge, phis_edge = jacobi_duals(ref_el, edge_weight, degree-4, degree) + Q_edge, phis_edge = jacobi_duals(ref_el, edge_weight, degree, degree-4) mapping = None self._reduced_dofs = {} @@ -75,7 +88,7 @@ def __init__(self, ref_el, degree): if dim == 1: Q_ref, Phis = Q_edge, phis_edge[:moment_degree+1] elif dim < sd: - Q_ref, Phis = self._reference_duals(ref_el, dim, degree, moment_degree) + Q_ref, Phis = dubiner_duals(ref_el, dim, degree, moment_degree) for entity in sorted(top[dim]): cur = len(nodes) @@ -131,18 +144,6 @@ def __init__(self, ref_el, degree): super().__init__(nodes, ref_el, entity_ids) - def _reference_duals(self, ref_el, dim, degree, moment_degree): - """Compute moments of trial space V against a subset of itself. - """ - if dim == 0: - return None, [] - - facet = ref_el.construct_subelement(dim) - V = ONPolynomialSet(facet, moment_degree, scale="orthonormal") - Q = create_quadrature(facet, degree + V.degree) - phis = V.tabulate(Q.get_points())[(0,)*dim] - return Q, phis - def _interior_duals(self, ref_el, degree): """Compute div-div and eps-eps moments of the trial space against an orthonormal bases for div(V_0) and eps(V_0). @@ -253,9 +254,12 @@ def __init__(self, ref_complex, degree): self._reduced_dofs = {} for dim in sorted(top): entity_ids[dim] = {} - if dim > 0 and dim < sd: - moment_degree = degree - (dim+1) - Q_ref, Phis = self._reference_duals(ref_complex, dim, degree, moment_degree) + + moment_degree = degree - (dim+1) + if dim == 1: + Q_ref, Phis = jacobi_duals(ref_complex, 0, degree, moment_degree) + elif dim > 0 and dim < sd: + Q_ref, Phis = dubiner_duals(ref_complex, dim, degree, moment_degree) self._reduced_dofs[dim] = None for entity in sorted(top[dim]): @@ -295,18 +299,104 @@ def __init__(self, ref_el, degree): super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") +class DivStokesDual(dual_set.DualSet): + def __init__(self, ref_el, degree): + self.degree = degree + nodes = [] + entity_ids = {} + top = ref_el.get_topology() + sd = ref_el.get_spatial_dimension() + Q = create_quadrature(ref_el, 2*degree) + + mapping = None + for dim in sorted(top): + entity_ids[dim] = {} + for entity in sorted(top[dim]): + entity_ids[dim][entity] = [] + + # Vertex dofs + for entity in top[0]: + pt, = ref_el.make_points(0, entity, degree) + nodes.append(PointEvaluation(ref_el, pt)) + + if sd == 3: + # Edge dofs + Q_ref, Phis = jacobi_duals(ref_el, 0, degree, degree-2) + for entity in top[1]: + Q_edge, phis = map_duals(ref_el, 1, entity, mapping, Q_ref, Phis) + nodes.extend(IntegralMoment(ref_el, Q_edge, phi) for phi in phis) + + # Interior dof + nodes.append(IntegralMoment(ref_el, Q, numpy.ones(Q.get_weights().shape))) + + # Throwaway dofs + B = Bernstein(ref_el, degree) + ids = B.entity_dofs() + + indices = [] + for dim in sorted(top): + for entity in sorted(top[dim]): + if sd == 3 and dim == 1: + continue + start = 1 if dim % sd == 0 else 0 + indices.extend(ids[dim][entity][start:]) + + interior = ids[sd][0][0] + phis = B.tabulate(0, Q.get_points())[(0,)*sd] + phis -= phis[[interior]] + nodes.extend(IntegralMoment(ref_el, Q, phi) for phi in phis[indices]) + + entity_ids[sd][0] = list(range(len(nodes))) + super().__init__(nodes, ref_el, entity_ids) + + def get_indices(self, restriction_domain, take_closure=True): + """Return the list of dofs with support on the given restriction domain. + Allows for reduced Demkowicz elements, excluding the exterior + derivative of the previous space in the de Rham complex. + + :arg restriction_domain: can be 'reduced', 'interior', 'vertex', + 'edge', 'face' or 'facet' + :kwarg take_closure: Are we taking the closure of the restriction domain? + """ + if restriction_domain == "reduced": + sd = self.ref_el.get_spatial_dimension() + top = self.ref_el.get_topology() + num_bfs = 1 + len(top[0]) + if sd == 3: + num_bfs += len(top[1]) * (self.degree-1) + return list(range(num_bfs)) + else: + return dual_set.DualSet.get_indices(self, restriction_domain, take_closure=take_closure) + + +class DivStokes(finite_element.CiarletElement): + def __init__(self, ref_el, degree): + sd = ref_el.get_spatial_dimension() + if degree < 2*sd-1: + raise ValueError(f"{type(self).__name__} elements only valid for k >= {2*sd-1}") + + poly_set = ONPolynomialSet(ref_el, degree, variant="bubble") + dual = DivStokesDual(ref_el, degree) + formdegree = sd # n-form + super().__init__(poly_set, dual, degree, formdegree) + + if __name__ == "__main__": import matplotlib.pyplot as plt import FIAT - dim = 2 - degree = 8 + dim = 3 + degree = 6 ref_el = FIAT.reference_element.symmetric_simplex(dim) - # fe = Stokes(ref_el, degree) - fe = MacroStokes(ref_el, degree) + fe = Stokes(ref_el, degree) + # fe = MacroStokes(ref_el, degree) Q = create_quadrature(fe.ref_complex, 2 * degree) Qpts, Qwts = Q.get_points(), Q.get_weights() + df = DivStokes(ref_el, degree-1) + df = FIAT.RestrictedElement(df, restriction_domain="reduced") + print(df.tabulate(0, ref_el.vertices)[(0,)*dim]) + family = type(fe).__name__ domains = (None, "reduced", "facet") fig, axes = plt.subplots(nrows=2, ncols=len(domains), figsize=(6*len(domains), 6*2))