diff --git a/FIAT/barycentric_interpolation.py b/FIAT/barycentric_interpolation.py index 25862bf2f..8b1fba4a4 100644 --- a/FIAT/barycentric_interpolation.py +++ b/FIAT/barycentric_interpolation.py @@ -30,7 +30,7 @@ def barycentric_interpolation(nodes, wts, dmat, pts, order=0): sp_simplify = numpy.vectorize(simplify) else: sp_simplify = lambda x: x - phi = numpy.add.outer(-nodes, pts.flatten()) + phi = numpy.add.outer(-nodes.flatten(), pts.flatten()) with numpy.errstate(divide='ignore', invalid='ignore'): numpy.reciprocal(phi, out=phi) numpy.multiply(phi, wts[:, None], out=phi) @@ -49,7 +49,7 @@ def barycentric_interpolation(nodes, wts, dmat, pts, order=0): def make_dmat(x): """Returns Lagrange differentiation matrix and barycentric weights associated with x[j].""" - dmat = numpy.add.outer(-x, x) + dmat = numpy.add.outer(-x.flatten(), x.flatten()) numpy.fill_diagonal(dmat, 1.0) wts = numpy.prod(dmat, axis=0) numpy.reciprocal(wts, out=wts) @@ -59,22 +59,38 @@ def make_dmat(x): class LagrangeLineExpansionSet(expansions.LineExpansionSet): - """Lagrange polynomial expansion set for given points the line.""" + """Lagrange polynomial expansion set for given points on the line.""" def __init__(self, ref_el, pts): + if ref_el.get_spatial_dimension() != 1: + raise Exception("Must have a line") + pts = numpy.asarray(pts) self.points = pts - self.x = numpy.array(pts, dtype="d").flatten() + self.cell_node_map = expansions.compute_cell_point_map(ref_el, pts, unique=False) self.dmats = [None for _ in self.cell_node_map] self.weights = [None for _ in self.cell_node_map] self.nodes = [None for _ in self.cell_node_map] + self.affine_mappings = {} for cell, ibfs in self.cell_node_map.items(): - self.nodes[cell] = self.x[ibfs] - self.dmats[cell], self.weights[cell] = make_dmat(self.nodes[cell]) + x = pts[ibfs] + if ref_el.is_trace(): + verts = ref_el.get_vertices_of_subcomplex(ref_el.topology[1][cell]) + A, = numpy.diff(verts, axis=0) + A /= numpy.linalg.norm(A) + b = -numpy.dot(numpy.sum(verts, axis=0)/2, A.T) + self.affine_mappings[cell] = (A, b) + x = numpy.add(numpy.dot(x, A.T), b) + self.nodes[cell] = x + self.dmats[cell], self.weights[cell] = make_dmat(x) self.degree = max(len(wts) for wts in self.weights)-1 self.recurrence_order = self.degree + 1 - super().__init__(ref_el) - self.continuity = None if len(self.x) == sum(len(xk) for xk in self.nodes) else "C0" + self.ref_el = ref_el + self.variant = None + self.scale = 1 + self.continuity = None if len(pts) == sum(len(xk) for xk in self.nodes) else "C0" + self._dmats_cache = {} + self._cell_node_map_cache = {} def get_num_members(self, n): return len(self.points) @@ -89,6 +105,12 @@ def get_dmats(self, degree, cell=0): return [self.dmats[cell].T] def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): + try: + A, b = self.affine_mappings[cell] + ref_pts = numpy.add(numpy.dot(pts.reshape(-1, A.shape[-1]), A.T), b) + pts = ref_pts.reshape(*pts.shape[:-1], -1) + except KeyError: + pass return barycentric_interpolation(self.nodes[cell], self.weights[cell], self.dmats[cell], pts, order=order) diff --git a/FIAT/bubble.py b/FIAT/bubble.py index 0de7b2566..6ecfeb826 100644 --- a/FIAT/bubble.py +++ b/FIAT/bubble.py @@ -14,8 +14,8 @@ class CodimBubble(RestrictedElement): """Bubbles of a certain codimension.""" - def __init__(self, ref_el, degree, codim): - element = Lagrange(ref_el, degree) + def __init__(self, ref_el, degree, codim, variant=None): + element = Lagrange(ref_el, degree, variant=variant) cell_dim = ref_el.get_dimension() assert cell_dim == max(element.entity_dofs().keys()) @@ -29,12 +29,12 @@ def __init__(self, ref_el, degree, codim): class Bubble(CodimBubble): """The bubble finite element: the dofs of the Lagrange FE in the interior of the cell""" - def __init__(self, ref_el, degree): - super().__init__(ref_el, degree, codim=0) + def __init__(self, ref_el, degree, variant=None): + super().__init__(ref_el, degree, codim=0, variant=variant) class FacetBubble(CodimBubble): """The facet bubble finite element: the dofs of the Lagrange FE in the interior of the facets""" - def __init__(self, ref_el, degree): - super().__init__(ref_el, degree, codim=1) + def __init__(self, ref_el, degree, variant=None): + super().__init__(ref_el, degree, codim=1, variant=variant) diff --git a/FIAT/discontinuous_lagrange.py b/FIAT/discontinuous_lagrange.py index fd08f8fa3..f806c9278 100644 --- a/FIAT/discontinuous_lagrange.py +++ b/FIAT/discontinuous_lagrange.py @@ -14,6 +14,7 @@ from FIAT.barycentric_interpolation import LagrangePolynomialSet, get_lagrange_points from FIAT.polynomial_set import mis from FIAT.check_format_variant import parse_lagrange_variant +from FIAT.hierarchical import Legendre def make_entity_permutations(dim, npoints): @@ -215,6 +216,8 @@ class DiscontinuousLagrange(finite_element.CiarletElement): macroelement for Scott-Vogelius. """ def __new__(cls, ref_el, degree, variant="equispaced"): + if variant and variant.startswith("integral"): + return Legendre(ref_el, degree, variant=variant) if degree == 0: splitting, _ = parse_lagrange_variant(variant, discontinuous=True) if splitting is None and not ref_el.is_macrocell(): diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index c09db974f..539c87553 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -277,15 +277,15 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations): parent_cell = ref_el.get_parent() if parent_cell is None: return nodes, ref_el, entity_ids, entity_permutations - parent_ids = {} + parent_top = parent_cell.get_topology() + parent_ids = {dim: {entity: [] for entity in parent_top[dim]} for dim in parent_top} parent_permutations = None parent_to_children = ref_el.get_parent_to_children() - if all(isinstance(node, functional.PointEvaluation) for node in nodes): + if all(isinstance(node, functional.PointEvaluation) for node in nodes) and not ref_el.is_trace(): # Merge Lagrange dual with lexicographical reordering parent_nodes = [] for dim in sorted(parent_to_children): - parent_ids[dim] = {} for entity in sorted(parent_to_children[dim]): cur = len(parent_nodes) for child_dim, child_entity in parent_to_children[dim][entity]: @@ -296,9 +296,7 @@ def merge_entities(nodes, ref_el, entity_ids, entity_permutations): # Merge everything else with the same node ordering parent_nodes = nodes for dim in sorted(parent_to_children): - parent_ids[dim] = {} for entity in sorted(parent_to_children[dim]): - parent_ids[dim][entity] = [] for child_dim, child_entity in parent_to_children[dim][entity]: parent_ids[dim][entity].extend(entity_ids[child_dim][child_entity]) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 69d0d4b5a..fd1de72b8 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -268,7 +268,7 @@ def __init__(self, ref_el, scale=None, variant=None): base_ref_el = reference_element.default_simplex(sd) base_verts = base_ref_el.get_vertices() - self.affine_mappings = [reference_element.make_affine_mapping( + self.affine_mappings = [get_affine_mapping( ref_el.get_vertices_of_subcomplex(top[sd][cell]), base_verts) for cell in top[sd]] if scale is None: @@ -354,6 +354,9 @@ def _tabulate(self, n, pts, order=0): if not self.ref_el.is_macrocell(): return phis[0] + if self.ref_el.is_trace(): + phis = trace_tabulation(self.ref_el, cell_point_map, order, pts, phis) + if pts.dtype == object: # If binning is undefined, scale by the characteristic function of each subcell Xi = compute_partition_of_unity(self.ref_el, pts, unique=unique) @@ -362,13 +365,14 @@ def _tabulate(self, n, pts, order=0): phi[alpha] *= Xi[cell] elif not unique: # If binning is not unique, divide by the multiplicity of each point - mult = numpy.zeros(pts.shape[:-1]) + mult = numpy.zeros(pts.shape[:-1], dtype=int) for cell, ipts in cell_point_map.items(): mult[ipts] += 1 - for cell, ipts in cell_point_map.items(): - phi = phis[cell] - for alpha in phi: - phi[alpha] /= mult[None, ipts] + if (mult != 1).any(): + for cell, ipts in cell_point_map.items(): + phi = phis[cell] + for alpha in phi: + phi[alpha] /= mult[None, ipts] # Insert subcell tabulations into the corresponding submatrices idx = lambda *args: args if args[1] is Ellipsis else numpy.ix_(*args) @@ -377,6 +381,9 @@ def _tabulate(self, n, pts, order=0): result = {} base_phi = tuple(phis.values())[0] for alpha in base_phi: + if isinstance(base_phi[alpha], Exception): + result[alpha] = base_phi[alpha] + continue dtype = base_phi[alpha].dtype result[alpha] = numpy.zeros((num_phis, *pts.shape[:-1]), dtype=dtype) for cell in cell_point_map: @@ -497,8 +504,7 @@ def get_dmats(self, degree, cell=0): def tabulate(self, n, pts): if len(pts) == 0: return numpy.array([]) - sd = self.ref_el.get_spatial_dimension() - return self._tabulate(n, pts)[(0,) * sd] + return tuple(self._tabulate(n, pts).values())[0] def tabulate_derivatives(self, n, pts): from FIAT.polynomial_set import mis @@ -591,13 +597,30 @@ def __init__(self, ref_el, **kwargs): super().__init__(ref_el, **kwargs) +def get_affine_mapping(xs, ys): + """Constructs (A,b) such that x --> A * x + b is the affine + mapping from the simplex defined by xs to the simplex defined by ys. + Uses least-squares when the simplex dimension does not match the spatial + dimension. + """ + X = numpy.asarray(xs) + Y = numpy.asarray(ys) + DX = X[1:] - X[:1] + DY = Y[1:] - Y[:1] + if DX.shape[0] == DX.shape[1]: + AT = numpy.linalg.solve(DX, DY) + else: + AT, *_ = numpy.linalg.lstsq(DX, DY) + b = Y[0] - numpy.dot(X[0], AT) + return AT.T, b + + def polynomial_dimension(ref_el, n, continuity=None): """Returns the dimension of the space of polynomials of degree no greater than n on the reference complex.""" if ref_el.get_shape() == reference_element.POINT: if n > 0: raise ValueError("Only degree zero polynomials supported on point elements.") - return 1 top = ref_el.get_topology() if continuity == "C0": space_dimension = sum(math.comb(n - 1, dim) * len(top[dim]) for dim in top) @@ -679,7 +702,10 @@ def compute_cell_point_map(ref_el, pts, unique=True, tol=1E-12): return {cell: Ellipsis for cell in sorted(top[sd])} # The distance to the nearest cell is equal to the distance to the parent cell - best = ref_el.get_parent().distance_to_point_l1(pts, rescale=True) + parent = ref_el + while parent.get_parent() is not None: + parent = parent.get_parent() + best = parent.distance_to_point_l1(pts, rescale=True) tol = best + tol cell_point_map = {} @@ -735,3 +761,29 @@ def compute_partition_of_unity(ref_el, pt, unique=True, tol=1E-12): mult = sum(masks) masks = [m / mult for m in masks] return masks + + +def trace_tabulation(ref_el, cell_point_map, order, pts, phis): + """Lift trace tabulations into the cells and raise TraceError on invalid tabulations.""" + from FIAT.polynomial_set import mis + from FIAT.hdiv_trace import TraceError + parent = ref_el.get_parent() + tdim = ref_el.get_spatial_dimension() + gdim = parent.get_spatial_dimension() + for cell in phis: + # Lift facet keys to cell keys + phi = phis[cell][(0,) * tdim] + # Raise TraceError on gradient tabulations + msg = "Gradients on trace elements are not well-defined." + phis[cell] = {alpha: phi if sum(alpha) == 0 else TraceError(msg) + for i in range(order+1) + for alpha in mis(gdim, i)} + + if sum(len(cell_point_map[cell]) for cell in cell_point_map) < len(pts): + # Raise TraceError when interior points fail to be binned on facets + for cell in parent.topology[gdim]: + msg = "The HDivTrace element can only be tabulated on facets." + phis[cell] = {alpha: TraceError(msg) + for i in range(order+1) + for alpha in mis(gdim, i)} + return phis diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index 264205a3e..3c26d0b47 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -4,21 +4,79 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later -import numpy as np -from FIAT.barycentric_interpolation import get_lagrange_points from FIAT.discontinuous_lagrange import DiscontinuousLagrange -from FIAT.dual_set import DualSet -from FIAT.finite_element import FiniteElement -from FIAT.functional import PointEvaluation -from FIAT.polynomial_set import mis -from FIAT.reference_element import (ufc_simplex, POINT, - LINE, QUADRILATERAL, - TRIANGLE, TETRAHEDRON, - TENSORPRODUCT) -from FIAT.tensor_product import TensorProductElement +from FIAT.reference_element import SimplicialComplex + + +class TraceSimplicialComplex(SimplicialComplex): + + def __init__(self, parent, codim=1): + sd = parent.get_spatial_dimension() - codim + parent_top = parent.topology + topology = dict(parent_top) + for dim in parent_top: + if dim > sd: + topology.pop(dim) + + facet = parent.construct_subelement(sd) + self.base_ref_el = facet + self._parent_complex = parent + self._parent_simplex = parent + + # dict mapping child facets to their parent facet + parent_to_children = {d: {e: [(d, e)] + for e in topology[d]} + for d in topology} + # dict mapping parent facets to their children + child_to_parent = {d: {e: (d, e) + for e in topology[d]} + for d in topology} + self._child_to_parent = child_to_parent + self._parent_to_children = parent_to_children + + # dict mapping cells to their boundary facets for each dimension, + # while respecting the ordering on the parent simplex + connectivity = {cell: {dim: [] for dim in topology} for cell in topology[sd]} + for cell in topology[sd]: + for dim in topology: + for entity in topology[dim]: + if set(topology[dim][entity]) <= set(topology[sd][cell]): + connectivity[cell][dim].append(entity) + self._cell_connectivity = connectivity + + # dict mapping subentity dimension to interior facets + interior_facets = {dim: [entity for entity in child_to_parent[dim] + if child_to_parent[dim][entity][0] == sd] + for dim in sorted(child_to_parent)} + self._interior_facets = interior_facets + super().__init__(facet.shape, parent.vertices, topology) + + def get_spatial_dimension(self): + return self.base_ref_el.get_spatial_dimension() + + def get_dimension(self): + return self.base_ref_el.get_dimension() + + def construct_subelement(self, dim): + return self._parent_simplex.construct_subelement(dim) + + def get_parent(self): + return self._parent_simplex + + def get_parent_complex(self): + return self._parent_complex + + def get_parent_to_children(self): + return self._parent_to_children + + def get_cell_connectivity(self): + return self._cell_connectivity + + def is_trace(self): + return True -# Numerical tolerance for facet-entity identifications -epsilon = 1e-10 + def is_macrocell(self): + return True class TraceError(Exception): @@ -30,7 +88,7 @@ def __init__(self, msg): self.msg = msg -class HDivTrace(FiniteElement): +class HDivTrace(DiscontinuousLagrange): """Class implementing the trace of hdiv elements. This class is a stand-alone element family that produces a DG-facet field. This element is what's produced after performing the trace @@ -39,358 +97,12 @@ class HDivTrace(FiniteElement): This element is also known as the discontinuous trace field that arises in several DG formulations. """ - - def __init__(self, ref_el, degree, variant=None): + def __new__(cls, ref_el, degree, variant="equispaced_interior"): """Constructor for the HDivTrace element. - :arg ref_el: A reference element, which may be a tensor product - cell. - :arg degree: The degree of approximation. If on a tensor product - cell, then provide a tuple of degrees if you want - varying degrees. + :arg ref_el: A reference element. + :arg degree: The degree of approximation. :arg variant: The point distribution variant passed on to recursivenodes. """ - sd = ref_el.get_spatial_dimension() - if sd in (0, 1): - raise ValueError("Cannot take the trace of a %d-dim cell." % sd) - - # Store the degrees if on a tensor product cell - if ref_el.get_shape() == TENSORPRODUCT: - try: - degree = tuple(degree) - except TypeError: - degree = (degree,) * len(ref_el.cells) - - assert len(ref_el.cells) == len(degree), ( - "Number of specified degrees must be equal to the number of cells." - ) - else: - if ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON, QUADRILATERAL]: - raise NotImplementedError( - "Trace element on a %s not implemented" % type(ref_el) - ) - # Cannot have varying degrees for these reference cells - if isinstance(degree, tuple): - raise ValueError("Must have a tensor product cell if providing multiple degrees") - - # Initialize entity dofs and construct the DG elements - # for the facets - facet_sd = sd - 1 - dg_elements = {} - entity_dofs = {} - topology = ref_el.get_topology() - for top_dim, entities in topology.items(): - cell = ref_el.construct_subelement(top_dim) - entity_dofs[top_dim] = {} - - # We have a facet entity! - if cell.get_spatial_dimension() == facet_sd: - dg_elements[top_dim] = construct_dg_element(cell, degree, variant) - # Initialize - for entity in entities: - entity_dofs[top_dim][entity] = [] - - # Compute the dof numbering for all facet entities - # and extract nodes - offset = 0 - pts = [] - for facet_dim in sorted(dg_elements): - element = dg_elements[facet_dim] - nf = element.space_dimension() - num_facets = len(topology[facet_dim]) - - facet_pts = get_lagrange_points(element.dual_basis()) - for i in range(num_facets): - entity_dofs[facet_dim][i] = list(range(offset, offset + nf)) - offset += nf - - # Run over nodes and collect the points for point evaluations - transform = ref_el.get_entity_transform(facet_dim, i) - pts.extend(transform(facet_pts)) - - # Setting up dual basis - only point evaluations - nodes = [PointEvaluation(ref_el, pt) for pt in pts] - dual = DualSet(nodes, ref_el, entity_dofs) - - # Degree of the element - deg = max([e.degree() for e in dg_elements.values()]) - - super().__init__(ref_el, dual, order=deg, - formdegree=facet_sd, - mapping="affine") - - # Set up facet elements - self.dg_elements = dg_elements - - # Degree for quadrature rule - self.polydegree = deg - - def degree(self): - """Return the degree of the (embedding) polynomial space.""" - return self.polydegree - - def get_nodal_basis(self): - """Return the nodal basis, encoded as a PolynomialSet object, - for the finite element.""" - raise NotImplementedError("get_nodal_basis not implemented for the trace element.") - - def get_coeffs(self): - """Return the expansion coefficients for the basis of the - finite element.""" - raise NotImplementedError("get_coeffs not implemented for the trace element.") - - def tabulate(self, order, points, entity=None): - """Return tabulated values of derivatives up to a given order of - basis functions at given points. - - :arg order: The maximum order of derivative. - :arg points: An iterable of points. - :arg entity: Optional (dimension, entity number) pair - indicating which topological entity of the - reference element to tabulate on. If ``None``, - tabulated values are computed by geometrically - approximating which facet the points are on. - - .. note :: - - Performing illegal tabulations on this element will result in either - a tabulation table of `numpy.nan` arrays (`entity=None` case), or - insertions of the `TraceError` exception class. This is due to the - fact that performing cell-wise tabulations, or asking for any order - of derivative evaluations, are not mathematically well-defined. - """ - sd = self.ref_el.get_spatial_dimension() - facet_sd = sd - 1 - - # Initializing dictionary with zeros - phivals = {} - for i in range(order + 1): - alphas = mis(sd, i) - for alpha in alphas: - phivals[alpha] = np.zeros(shape=(self.space_dimension(), len(points))) - - evalkey = (0,) * sd - - # If entity is None, identify facet using numerical tolerance and - # return the tabulated values - if entity is None: - # NOTE: Numerical approximation of the facet id is currently only - # implemented for simplex reference cells. - if self.ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON]: - raise NotImplementedError( - "Tabulating this element on a %s cell without providing " - "an entity is not currently supported." % type(self.ref_el) - ) - - # Attempt to identify which facet (if any) the given points are on - vertices = self.ref_el.vertices - coordinates = barycentric_coordinates(points, vertices) - unique_facet, success = extract_unique_facet(coordinates) - - # If not successful, return NaNs - if not success: - for key in phivals: - phivals[key] = np.full(shape=(self.space_dimension(), len(points)), fill_value=np.nan) - - return phivals - - # Otherwise, extract non-zero values and insertion indices - else: - # Map points to the reference facet - new_points = map_to_reference_facet(points, vertices, unique_facet) - - # Retrieve values by tabulating the DG element - element = self.dg_elements[facet_sd] - nf = element.space_dimension() - nonzerovals, = element.tabulate(order, new_points).values() - indices = slice(nf * unique_facet, nf * (unique_facet + 1)) - - else: - entity_dim, _ = entity - - # If the user is directly specifying cell-wise tabulation, return - # TraceErrors in dict for appropriate handling in the form compiler - if entity_dim not in self.dg_elements: - for key in phivals: - msg = "The HDivTrace element can only be tabulated on facets." - phivals[key] = TraceError(msg) - - return phivals - - else: - # Retrieve function evaluations (order = 0 case) - offset = 0 - for facet_dim in sorted(self.dg_elements): - element = self.dg_elements[facet_dim] - nf = element.space_dimension() - num_facets = len(self.ref_el.get_topology()[facet_dim]) - - # Loop over the number of facets until we find a facet - # with matching dimension and id - for i in range(num_facets): - # Found it! Grab insertion indices - if (facet_dim, i) == entity: - nonzerovals, = element.tabulate(0, points).values() - indices = slice(offset, offset + nf) - - offset += nf - - # If asking for gradient evaluations, insert TraceError in - # gradient slots - if order > 0: - msg = "Gradients on trace elements are not well-defined." - for key in phivals: - if key != evalkey: - phivals[key] = TraceError(msg) - - # Insert non-zero values in appropriate place - phivals[evalkey][indices, :] = nonzerovals - - return phivals - - def value_shape(self): - """Return the value shape of the finite element functions.""" - return () - - def dmats(self): - """Return dmats: expansion coefficients for basis function - derivatives.""" - raise NotImplementedError("dmats not implemented for the trace element.") - - def get_num_members(self, arg): - """Return number of members of the expansion set.""" - raise NotImplementedError("get_num_members not implemented for the trace element.") - - @staticmethod - def is_nodal(): - return True - - -def construct_dg_element(ref_el, degree, variant): - """Constructs a discontinuous galerkin element of a given degree - on a particular reference cell. - """ - if ref_el.get_shape() in [LINE, TRIANGLE]: - dg_element = DiscontinuousLagrange(ref_el, degree, variant) - - # Quadrilateral facets could be on a FiredrakeQuadrilateral. - # In this case, we treat this as an interval x interval cell: - elif ref_el.get_shape() == QUADRILATERAL: - dg_line = DiscontinuousLagrange(ufc_simplex(1), degree, variant) - dg_element = TensorProductElement(dg_line, dg_line) - - # This handles the more general case for facets: - elif ref_el.get_shape() == TENSORPRODUCT: - assert len(degree) == len(ref_el.cells), ( - "Must provide the same number of degrees as the number " - "of cells that make up the tensor product cell." - ) - sub_elements = [construct_dg_element(c, d, variant) - for c, d in zip(ref_el.cells, degree) - if c.get_shape() != POINT] - - if len(sub_elements) > 1: - dg_element = TensorProductElement(*sub_elements) - else: - dg_element, = sub_elements - - else: - raise NotImplementedError( - "Reference cells of type %s not currently supported" % type(ref_el) - ) - - return dg_element - - -# The following functions are credited to Marie E. Rognes: -def extract_unique_facet(coordinates, tolerance=epsilon): - """Determines whether a set of points (described in its barycentric coordinates) - are all on one of the facet sub-entities, and return the particular facet and - whether the search has been successful. - - :arg coordinates: A set of points described in barycentric coordinates. - :arg tolerance: A fixed tolerance for geometric identifications. - """ - facets = [] - for c in coordinates: - on_facet = set([i for (i, l) in enumerate(c) if abs(l) < tolerance]) - facets += [on_facet] - - unique_facet = facets[0] - for f in facets: - unique_facet = unique_facet & f - - # Handle coordinates not on facets - if len(unique_facet) != 1: - return (None, False) - - # If we have a unique facet, return it and success - return (unique_facet.pop(), True) - - -def barycentric_coordinates(points, vertices): - """Computes the barycentric coordinates for a set of points relative to a - simplex defined by a set of vertices. - - :arg points: A set of points. - :arg vertices: A set of vertices that define the simplex. - """ - - # Form mapping matrix - T = (np.asarray(vertices[:-1]) - vertices[-1]).T - invT = np.linalg.inv(T) - - points = np.asarray(points) - bary = np.einsum("ij,kj->ki", invT, (points - vertices[-1])) - last = (1 - bary.sum(axis=1)) - return np.concatenate([bary, last[..., np.newaxis]], axis=1) - - -def map_from_reference_facet(point, vertices): - """Evaluates the physical coordinate of a point using barycentric - coordinates. - - :arg point: The reference points to be mapped to the facet. - :arg vertices: The vertices defining the physical element. - """ - - # Compute the barycentric coordinates of the point relative to the reference facet - reference_simplex = ufc_simplex(len(vertices) - 1) - reference_vertices = reference_simplex.get_vertices() - coords = barycentric_coordinates([point, ], reference_vertices)[0] - - # Evaluates the physical coordinate of the point using barycentric coordinates - point = sum(vertices[j] * coords[j] for j in range(len(coords))) - return tuple(point) - - -def map_to_reference_facet(points, vertices, facet): - """Given a set of points and vertices describing a facet of a simplex in n-dimensional - coordinates (where the points lie on the facet), map the points to the reference simplex - of dimension (n-1). - - :arg points: A set of points in n-D. - :arg vertices: A set of vertices describing a facet of a simplex in n-D. - :arg facet: Integer representing the facet number. - """ - - # Compute the barycentric coordinates of the points with respect to the - # full physical simplex - all_coords = barycentric_coordinates(points, vertices) - - # Extract vertices of the reference facet - reference_facet_simplex = ufc_simplex(len(vertices) - 2) - reference_vertices = reference_facet_simplex.get_vertices() - - reference_points = [] - for (i, coords) in enumerate(all_coords): - # Extract the correct subset of barycentric coordinates since we know - # which facet we are on - new_coords = [coords[j] for j in range(len(coords)) if j != facet] - - # Evaluate the reference coordinate of a point in barycentric coordinates - reference_pt = sum(np.asarray(reference_vertices[j]) * new_coords[j] - for j in range(len(new_coords))) - - reference_points += [reference_pt] - return reference_points + facets = TraceSimplicialComplex(ref_el) + return DiscontinuousLagrange(facets, degree, variant=variant) diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py index 17c415710..535ea7d7b 100644 --- a/FIAT/hierarchical.py +++ b/FIAT/hierarchical.py @@ -51,13 +51,14 @@ def __init__(self, ref_el, degree, codim=0): continue ref_facet = ref_el.construct_subelement(dim) - poly_set = ONPolynomialSet(ref_facet, degree) + poly_set = ONPolynomialSet(ref_facet, degree, scale=1) Q_ref = create_quadrature(ref_facet, 2 * degree) phis = poly_set.tabulate(Q_ref.get_points())[(0,) * dim] for entity in sorted(top[dim]): cur = len(nodes) + Jdet = ref_el.volume_of_subcomplex(dim, entity) Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref) - nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi) for phi in phis) + nodes.extend(functional.IntegralMoment(ref_el, Q_facet, phi/Jdet) for phi in phis) entity_ids[dim][entity] = list(range(cur, len(nodes))) entity_permutations[dim][entity] = perms @@ -69,7 +70,7 @@ class Legendre(finite_element.CiarletElement): def __new__(cls, ref_el, degree, variant=None): if degree == 0: splitting, _ = parse_lagrange_variant(variant, integral=True) - if splitting is None: + if splitting is None and not ref_el.is_macrocell(): # FIXME P0 on the split requires implementing SplitSimplicialComplex.symmetry_group_size() return P0.P0(ref_el) return super().__new__(cls) diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index eead0cfa6..d8219a4ce 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -10,6 +10,7 @@ from FIAT.barycentric_interpolation import LagrangePolynomialSet, get_lagrange_points from FIAT.reference_element import LINE from FIAT.check_format_variant import parse_lagrange_variant +from FIAT.hierarchical import IntegratedLegendre class LagrangeDualSet(dual_set.DualSet): @@ -72,6 +73,11 @@ class Lagrange(finite_element.CiarletElement): entity id. The DOFs are always sorted by the entity ordering and then lexicographically by lattice multiindex. """ + def __new__(cls, ref_el, degree, variant="equispaced", sort_entities=False): + if variant and variant.startswith("integral"): + return IntegratedLegendre(ref_el, degree, variant=variant) + return super().__new__(cls) + def __init__(self, ref_el, degree, variant="equispaced", sort_entities=False): splitting, point_variant = parse_lagrange_variant(variant) if splitting is not None: diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 534c18642..2d158fce4 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -68,8 +68,11 @@ def tabulate_new(self, pts): def tabulate(self, pts, jet_order=0): """Returns the values of the polynomial set.""" - base_vals = self.expansion_set._tabulate(self.embedded_degree, pts, order=jet_order) - result = {alpha: numpy.dot(self.coeffs, base_vals[alpha]) for alpha in base_vals} + result = self.expansion_set._tabulate(self.embedded_degree, pts, order=jet_order) + for alpha in result: + if isinstance(result[alpha], Exception): + continue + result[alpha] = numpy.dot(self.coeffs, result[alpha]) return result def get_expansion_set(self): diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 547da2c04..553658acd 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -94,6 +94,8 @@ def make_lattice(verts, n, interior=0, variant=None): family = _decode_family(family) D = len(verts) X = numpy.array(verts) + if D == 1 and n == 0: + return list(map(tuple, X)) get_point = lambda alpha: tuple(numpy.dot(_recursive(D - 1, n, alpha, family), X)) return list(map(get_point, multiindex_equal(D, n, interior))) @@ -323,6 +325,9 @@ def extrinsic_orientation_permutation_map(self): def is_simplex(self): return False + def is_trace(self): + return False + def is_macrocell(self): return False @@ -341,6 +346,10 @@ def get_parent_complex(self): def is_parent(self, other, strict=False): """Return whether this cell is the parent of the other cell.""" parent = other + while self.is_trace(): + self = self.get_parent_complex() + while parent.is_trace(): + parent = parent.get_parent_complex() if strict: parent = parent.get_parent_complex() while parent is not None: @@ -357,6 +366,8 @@ def __eq__(self, other): return False atop = self.get_topology() btop = other.get_topology() + if max(atop) != max(btop): + return False for dim in atop: if set(atop[dim].values()) != set(btop[dim].values()): return False @@ -566,7 +577,6 @@ def get_entity_transform(self, dim, entity): """ topology = self.get_topology() celldim = self.get_spatial_dimension() - codim = celldim - dim if dim == 0: # Special case vertices. i, = topology[dim][entity] @@ -578,7 +588,7 @@ def get_entity_transform(self, dim, entity): else: subcell = self.construct_subelement(dim) subdim = subcell.get_spatial_dimension() - assert subdim == celldim - codim + assert subdim == dim # Entity vertices in entity space. v_e = numpy.asarray(subcell.get_vertices()) @@ -608,30 +618,34 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): if len(points) == 0: return points if entity is None: - entity = (self.get_spatial_dimension(), 0) + entity = (self.get_dimension(), 0) entity_dim, entity_id = entity - top = self.get_topology() - sd = self.get_spatial_dimension() - - # get a subcell containing the entity and the restriction indices of the entity - indices = slice(None) - subcomplex = top[entity_dim][entity_id] - if entity_dim != sd: - cell_id = self.connectivity[(entity_dim, sd)][0][0] - indices = [i for i, v in enumerate(top[sd][cell_id]) if v in subcomplex] - subcomplex = top[sd][cell_id] + subcomplex = self.topology[entity_dim][entity_id] + + tdim = sum(entity_dim) if isinstance(entity_dim, tuple) else entity_dim + gdim = len(self.vertices[0]) + if tdim != gdim: + # Add extra vertices to construct a subcell containing the entity + parent = self.get_parent_complex() if self.is_trace() else self + cell_dim = parent.get_dimension() + cell_id = min(parent.connectivity[(entity_dim, cell_dim)][entity_id]) + subcell = parent.topology[cell_dim][cell_id] + subcomplex += tuple(set(subcell) - set(subcomplex))[:gdim-tdim] cell_verts = self.get_vertices_of_subcomplex(subcomplex) - ref_verts = numpy.eye(sd + 1) + ref_verts = numpy.eye(len(cell_verts)) A, b = make_affine_mapping(cell_verts, ref_verts) - A, b = A[indices], b[indices] if rescale: - # rescale barycentric coordinates by the height wrt. to the facet + # rescale barycentric coordinates by the height w.r.t. the facet h = 1 / numpy.linalg.norm(A, axis=1) b *= h A *= h[:, None] - out = numpy.dot(points, A.T) - return numpy.add(out, b, out=out) + bary = numpy.dot(points, A.T) + numpy.add(bary, b, out=bary) + if tdim != gdim: + # the extra coordinates are always negated to mark points outside the facet + bary[:, tdim-gdim:] = -abs(bary[:, tdim-gdim:]) + return bary def compute_bubble(self, points, entity=None): """Returns the lowest-order bubble on an entity evaluated at the given @@ -1249,9 +1263,9 @@ def distance_to_point_l1(self, point, rescale=False): For more information see the docstring for the UFCSimplex method.""" subcell_dimensions = self.get_dimension() - assert len(point) == sum(subcell_dimensions) - point_slices = TensorProductCell._split_slices(subcell_dimensions) point = numpy.asarray(point) + assert point.shape[-1] == sum(subcell_dimensions) + point_slices = TensorProductCell._split_slices(subcell_dimensions) return sum(c.distance_to_point_l1(point[..., s], rescale=rescale) for c, s in zip(self.cells, point_slices)) @@ -1390,6 +1404,17 @@ def construct_subelement(self, dimension): sub_element = self.product.construct_subelement((dimension,) + (0,)*(len(self.product.cells) - 1)) return flatten_reference_cube(sub_element) + def get_facet_element(self): + sd = self.get_spatial_dimension() + return self.construct_subelement(sd-1) + + def make_points(self, dim, entity, order, variant=None, interior=1): + if dim == 1: + entity_verts = self.get_vertices_of_subcomplex(self.topology[dim][entity]) + return make_lattice(entity_verts, order, interior=interior, variant=variant) + else: + raise NotImplementedError("Cannot make points on non-simplex facets yet.") + def get_entity_transform(self, dim, entity_i): """Returns a mapping of point coordinates from the `entity_i`-th subentity of dimension `dim` to the cell. diff --git a/finat/__init__.py b/finat/__init__.py index 5390fa779..dde090cb0 100644 --- a/finat/__init__.py +++ b/finat/__init__.py @@ -32,7 +32,7 @@ from .cube import FlattenedDimensions # noqa: F401 from .discontinuous import DiscontinuousElement # noqa: F401 from .enriched import EnrichedElement # noqa: F401 -from .hdivcurl import HCurlElement, HDivElement # noqa: F401 +from .hdivcurl import HCurlElement, HDivElement, HDivTraceElement # noqa: F401 from .mixed import MixedElement # noqa: F401 from .nodal_enriched import NodalEnrichedElement # noqa: F401 from .quadrature_element import QuadratureElement, make_quadrature_element # noqa: F401 diff --git a/finat/element_factory.py b/finat/element_factory.py index 8e038dccd..bb665849d 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -158,19 +158,21 @@ def convert_finiteelement(element, **kwargs): codim = 1 if element.family() == "Boundary Quadrature" else 0 return finat.make_quadrature_element(cell, degree, scheme, codim), set() lmbda = supported_elements[element.family()] - if element.family() == "Real" and element.cell.cellname() in {"quadrilateral", "hexahedron"}: - lmbda = None - element = finat.ufl.FiniteElement("DQ", element.cell, 0) + if element.cell.cellname() in supported_tensor_product_cells: + if element.family() == "Real": + lmbda = None + element = finat.ufl.FiniteElement("DQ", element.cell, 0) + elif element.family() == "HDiv Trace": + lmbda = None + if lmbda is None: - if element.cell.cellname() == "quadrilateral": - # Handle quadrilateral short names like RTCF and RTCE. - element = element.reconstruct(cell=quadrilateral_tpc) - elif element.cell.cellname() == "hexahedron": - # Handle hexahedron short names like NCF and NCE. - element = element.reconstruct(cell=hexahedron_tpc) - else: + # Handle quadrilateral short names like RTCE/F and NCE/F. + try: + tpc = supported_tensor_product_cells[element.cell.cellname()] + except KeyError: raise ValueError("%s is supported, but handled incorrectly" % element.family()) + element = element.reconstruct(cell=tpc) finat_elem, deps = _create_element(element, **kwargs) return finat.FlattenedDimensions(finat_elem), deps @@ -219,7 +221,7 @@ def convert_finiteelement(element, **kwargs): lmbda = finat.DiscontinuousLagrange finat_kwargs["variant"] = kind - elif element.family() == "HDiv Trace": + elif element.family() in {"HDiv Trace", "Bubble", "FacetBubble"}: finat_kwargs["variant"] = kind elif element.variant() is not None: @@ -285,6 +287,12 @@ def convert_tensorproductelement(element, **kwargs): return finat.TensorProductElement(elements), deps +@convert.register(finat.ufl.HDivTraceElement) +def convert_hdivtraceelement(element, **kwargs): + finat_elem, deps = _create_element(element._element, **kwargs) + return finat.HDivTraceElement(finat_elem), deps + + @convert.register(finat.ufl.HDivElement) def convert_hdivelement(element, **kwargs): finat_elem, deps = _create_element(element._element, **kwargs) @@ -310,6 +318,11 @@ def convert_restrictedelement(element, **kwargs): hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval) quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval) +supported_tensor_product_cells = { + "quadrilateral": quadrilateral_tpc, + "hexahedron": hexahedron_tpc, +} + _cache = weakref.WeakKeyDictionary() diff --git a/finat/fiat_elements.py b/finat/fiat_elements.py index c4799b1ae..b94cc1a7e 100644 --- a/finat/fiat_elements.py +++ b/finat/fiat_elements.py @@ -112,7 +112,8 @@ def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None): index_shape = (self._element.space_dimension(),) for alpha, fiat_table in fiat_result.items(): if isinstance(fiat_table, Exception): - result[alpha] = gem.Failure(self.index_shape + self.value_shape, fiat_table) + shape = ps.points.shape[:-1] + self.index_shape + self.value_shape + result[alpha] = gem.partial_indexed(gem.Failure(shape, fiat_table), ps.indices) continue derivative = sum(alpha) @@ -348,13 +349,13 @@ def __init__(self, cell, degree): class Bubble(ScalarFiatElement): - def __init__(self, cell, degree): - super().__init__(FIAT.Bubble(cell, degree)) + def __init__(self, cell, degree, variant=None): + super().__init__(FIAT.Bubble(cell, degree, variant=variant)) class FacetBubble(ScalarFiatElement): - def __init__(self, cell, degree): - super().__init__(FIAT.FacetBubble(cell, degree)) + def __init__(self, cell, degree, variant=None): + super().__init__(FIAT.FacetBubble(cell, degree, variant=variant)) class CrouzeixRaviart(ScalarFiatElement): diff --git a/finat/hdivcurl.py b/finat/hdivcurl.py index 3ca697744..c0bdb4454 100644 --- a/finat/hdivcurl.py +++ b/finat/hdivcurl.py @@ -2,6 +2,7 @@ from FIAT.reference_element import LINE import gem +import numpy from gem.utils import cached_property from finat.finiteelementbase import FiniteElementBase from finat.tensor_product import TensorProductElement @@ -58,7 +59,7 @@ def index_shape(self): def value_shape(self): return (self.cell.get_spatial_dimension(),) - def _transform_evaluation(self, core_eval): + def _transform_evaluation(self, core_eval, entity=None): beta = self.get_indices() zeta = self.get_value_indices() @@ -72,11 +73,11 @@ def promote(table): def basis_evaluation(self, order, ps, entity=None, coordinate_mapping=None): core_eval = self.wrappee.basis_evaluation(order, ps, entity) - return self._transform_evaluation(core_eval) + return self._transform_evaluation(core_eval, entity) def point_evaluation(self, order, refcoords, entity=None, coordinate_mapping=None): core_eval = self.wrappee.point_evaluation(order, refcoords, entity) - return self._transform_evaluation(core_eval) + return self._transform_evaluation(core_eval, entity) @property def dual_basis(self): @@ -92,6 +93,59 @@ def dual_basis(self): return gem.ComponentTensor(Q[zeta], beta + zeta), x +class HDivTraceElement(WrapperElementBase): + """HDiv Trace wrapper element for tensor product elements.""" + + def __init__(self, wrappee): + assert isinstance(wrappee, TensorProductElement) + if any(fe.formdegree is None for fe in wrappee.factors): + raise ValueError("Form degree of subelement is None, cannot HDiv Trace!") + + formdegree = sum(fe.formdegree for fe in wrappee.factors) + if formdegree != wrappee.cell.get_spatial_dimension() - 1: + raise ValueError("HDiv Trace requires (n-1)-form element!") + + self.support_dim = tuple(fe.formdegree for fe in wrappee.factors) + + super().__init__(wrappee, None) + + def _transform_evaluation(self, core_eval, entity): + if entity is None: + entity = (self.cell.get_dimension(), 0) + entity_dim, entity_id = entity + + if entity_dim == self.support_dim or sum(entity_dim) != sum(self.support_dim): + return core_eval + + def zero_failure(expr): + if isinstance(expr, gem.Failure): + return gem.Literal(numpy.zeros(expr.shape, expr.dtype)) + return expr.reconstruct(*map(zero_failure, expr.children)) + + # Create a zero tabulation with the same tensor-product structure + return {alpha: zero_failure(table) for alpha, table in core_eval.items()} + + @property + def formdegree(self): + return self.wrappee.formdegree + + @property + def value_shape(self): + return self.wrappee.value_shape + + @cached_property + def fiat_equivalent(self): + return self.wrappee.fiat_equivalent + + @property + def mapping(self): + return self.wrappee.mapping + + @property + def dual_basis(self): + return self.wrappee.dual_basis + + class HDivElement(WrapperElementBase): """H(div) wrapper element for tensor product elements.""" diff --git a/finat/ufl/__init__.py b/finat/ufl/__init__.py index 21a7d13d1..f36ca39db 100644 --- a/finat/ufl/__init__.py +++ b/finat/ufl/__init__.py @@ -15,7 +15,7 @@ from finat.ufl.enrichedelement import EnrichedElement, NodalEnrichedElement # noqa: F401 from finat.ufl.finiteelement import FiniteElement # noqa: F401 from finat.ufl.finiteelementbase import FiniteElementBase # noqa: F401 -from finat.ufl.hdivcurl import HCurlElement, HDivElement, WithMapping, HDiv, HCurl # noqa: F401 +from finat.ufl.hdivcurl import HCurlElement, HDivElement, HDivTraceElement, WithMapping, HDiv, HCurl # noqa: F401 from finat.ufl.mixedelement import MixedElement, TensorElement, VectorElement # noqa: F401 from finat.ufl.restrictedelement import RestrictedElement # noqa: F401 from finat.ufl.tensorproductelement import TensorProductElement # noqa: F401 diff --git a/finat/ufl/enrichedelement.py b/finat/ufl/enrichedelement.py index a054885bb..bf686c694 100644 --- a/finat/ufl/enrichedelement.py +++ b/finat/ufl/enrichedelement.py @@ -17,7 +17,7 @@ class EnrichedElementBase(FiniteElementBase): """The vector sum of several finite element spaces.""" - def __init__(self, *elements): + def __init__(self, *elements, family=None): """Doc.""" self._elements = elements @@ -49,9 +49,11 @@ def __init__(self, *elements): # Get name of subclass: EnrichedElement or NodalEnrichedElement class_name = self.__class__.__name__ + if family is None: + family = class_name # Initialize element data - FiniteElementBase.__init__(self, class_name, cell, degree, + FiniteElementBase.__init__(self, family, cell, degree, quad_scheme, reference_value_shape) def mapping(self): @@ -87,7 +89,8 @@ def variant(self): def reconstruct(self, **kwargs): """Doc.""" - return type(self)(*[e.reconstruct(**kwargs) for e in self._elements]) + family = kwargs.pop("family", self.family()) + return type(self)(*(e.reconstruct(**kwargs) for e in self._elements), family=family) @property def embedded_subdegree(self): diff --git a/finat/ufl/finiteelement.py b/finat/ufl/finiteelement.py index 4cf34c90d..0968f7e18 100644 --- a/finat/ufl/finiteelement.py +++ b/finat/ufl/finiteelement.py @@ -39,6 +39,7 @@ def __new__(cls, from finat.ufl.enrichedelement import EnrichedElement from finat.ufl.hdivcurl import HCurlElement as HCurl from finat.ufl.hdivcurl import HDivElement as HDiv + from finat.ufl.hdivcurl import HDivTraceElement as HDivTrace from finat.ufl.tensorproductelement import TensorProductElement family, short_name, degree, reference_value_shape, sobolev_space, mapping, embedded_degree = \ @@ -94,6 +95,24 @@ def __new__(cls, return EnrichedElement(HCurl(TensorProductElement(Qc_elt, Id_elt, cell=cell)), HCurl(TensorProductElement(Qd_elt, Ic_elt, cell=cell))) + elif family == "HDiv Trace": + cell_h, cell_v = cell.sub_cells() + if not isinstance(degree, tuple): + degree = (degree, degree) + hdegree, vdegree = degree + + dg_family = lambda cell: "DG" if cell.is_simplex() else "DQ" + is_interval = lambda cell: cell.cellname() == "interval" + + dg_h = FiniteElement(dg_family(cell_h), cell_h, hdegree, variant=variant) + tr_h = FiniteElement("HDiv Trace", cell_h, 0 if is_interval(cell_h) else hdegree, variant=variant) + + dg_v = FiniteElement(dg_family(cell_v), cell_v, vdegree, variant=variant) + tr_v = FiniteElement("HDiv Trace", cell_v, 0 if is_interval(cell_v) else vdegree, variant=variant) + + return EnrichedElement(HDivTrace(TensorProductElement(tr_h, dg_v, cell=cell)), + HDivTrace(TensorProductElement(dg_h, tr_v, cell=cell)), family=family) + elif family == "Q": return TensorProductElement(*[FiniteElement("CG", c, degree, variant=variant) for c in cell.sub_cells()], diff --git a/finat/ufl/hdivcurl.py b/finat/ufl/hdivcurl.py index ed188038e..0a80e397d 100644 --- a/finat/ufl/hdivcurl.py +++ b/finat/ufl/hdivcurl.py @@ -35,40 +35,38 @@ def __call__(self, element): HDiv = CallableSobolevSpace(HDivSobolevSpace.name, HDivSobolevSpace.parents) -class HDivElement(FiniteElementBase): - """A div-conforming version of an outer product element, assuming this makes mathematical sense.""" +class WrapperElementBase(FiniteElementBase): + """A modified version of a tensor product element.""" __slots__ = ("_element", ) - def __init__(self, element): - """Doc.""" + def __init__(self, element, reference_value_shape, sobolev_space, mapping): self._element = element + self._sobolev_space = sobolev_space + self._mapping = mapping - family = "TensorProductElement" + family = element.family() cell = element.cell degree = element.degree() quad_scheme = element.quadrature_scheme() - reference_value_shape = (element.cell.topological_dimension(),) - # Skipping TensorProductElement constructor! Bad code smell, refactor to avoid this non-inheritance somehow. - FiniteElementBase.__init__(self, family, cell, degree, - quad_scheme, reference_value_shape) + super().__init__(family, cell, degree, quad_scheme, reference_value_shape) def __repr__(self): """Doc.""" - return f"HDivElement({repr(self._element)})" + return f"{type(self).__name__}({repr(self._element)})" def mapping(self): """Doc.""" - return "contravariant Piola" + return self._mapping @property def sobolev_space(self): """Return the underlying Sobolev space.""" - return HDivSobolevSpace + return self._sobolev_space def reconstruct(self, **kwargs): """Doc.""" - return HDivElement(self._element.reconstruct(**kwargs)) + return type(self)(self._element.reconstruct(**kwargs)) def variant(self): """Doc.""" @@ -76,11 +74,11 @@ def variant(self): def __str__(self): """Doc.""" - return f"HDivElement({repr(self._element)})" + return f"{type(self).__name__}({repr(self._element)})" def shortstr(self): """Format as string for pretty printing.""" - return f"HDivElement({self._element.shortstr()})" + return f"HDivTraceElement({self._element.shortstr()})" @property def embedded_subdegree(self): @@ -93,63 +91,32 @@ def embedded_superdegree(self): return self._element.embedded_superdegree -class HCurlElement(FiniteElementBase): - """A curl-conforming version of an outer product element, assuming this makes mathematical sense.""" - __slots__ = ("_element",) +class HDivTraceElement(WrapperElementBase): + """A HDiv Trace version of a tensor product element, assuming this makes mathematical sense.""" + __slots__ = ("_element", ) def __init__(self, element): - """Doc.""" - self._element = element + reference_value_shape = () + super().__init__(element, reference_value_shape, L2, "identity") - family = "TensorProductElement" - cell = element.cell - degree = element.degree() - quad_scheme = element.quadrature_scheme() - cell = element.cell - reference_value_shape = (cell.topological_dimension(),) # TODO: Is this right? - # Skipping TensorProductElement constructor! Bad code smell, - # refactor to avoid this non-inheritance somehow. - FiniteElementBase.__init__(self, family, cell, degree, quad_scheme, - reference_value_shape) - - def __repr__(self): - """Doc.""" - return f"HCurlElement({repr(self._element)})" - def mapping(self): - """Doc.""" - return "covariant Piola" +class HDivElement(WrapperElementBase): + """A div-conforming version of a tensor product element, assuming this makes mathematical sense.""" + __slots__ = ("_element", ) - @property - def sobolev_space(self): - """Return the underlying Sobolev space.""" - return HCurlSobolevSpace + def __init__(self, element): + reference_value_shape = (element.cell.topological_dimension(),) + super().__init__(element, reference_value_shape, HDiv, "contravariant Piola") - def reconstruct(self, **kwargs): - """Doc.""" - return HCurlElement(self._element.reconstruct(**kwargs)) - def variant(self): - """Doc.""" - return self._element.variant() +class HCurlElement(WrapperElementBase): + """A curl-conforming version of a tensor product element, assuming this makes mathematical sense.""" + __slots__ = ("_element",) - def __str__(self): + def __init__(self, element): """Doc.""" - return f"HCurlElement({repr(self._element)})" - - def shortstr(self): - """Format as string for pretty printing.""" - return f"HCurlElement({self._element.shortstr()})" - - @property - def embedded_subdegree(self): - """Return embedded subdegree.""" - return self._element.embedded_subdegree - - @property - def embedded_superdegree(self): - """Return embedded superdegree.""" - return self._element.embedded_superdegree + reference_value_shape = (element.cell.topological_dimension(),) + super().__init__(element, reference_value_shape, HCurl, "covariant Piola") class WithMapping(FiniteElementBase): diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py index ffb192afc..d2586f57e 100644 --- a/test/FIAT/unit/test_fiat.py +++ b/test/FIAT/unit/test_fiat.py @@ -34,7 +34,7 @@ from FIAT.nedelec import Nedelec # noqa: F401 from FIAT.nedelec_second_kind import NedelecSecondKind # noqa: F401 from FIAT.regge import Regge # noqa: F401 -from FIAT.hdiv_trace import HDivTrace, map_to_reference_facet # noqa: F401 +from FIAT.hdiv_trace import HDivTrace # noqa: F401 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson # noqa: F401 from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlFirstKind # noqa: F401 from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlSecondKind # noqa: F401 @@ -59,7 +59,7 @@ from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1 # noqa: F401 from FIAT.guzman_neilan import GuzmanNeilanSecondKindH1 # noqa: F401 from FIAT.johnson_mercier import JohnsonMercier # noqa: F401 -from FIAT.bubble import Bubble +from FIAT.bubble import Bubble, FacetBubble # noqa: F401 from FIAT.enriched import EnrichedElement # noqa: F401 from FIAT.nodal_enriched import NodalEnrichedElement from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen # noqa: F401 @@ -320,8 +320,15 @@ def __init__(self, a, b): "Histopolation(I, 1)", "Histopolation(I, 2)", "Bubble(I, 2)", + "Bubble(I, 2, 'integral')", "Bubble(T, 3)", + "Bubble(T, 3, 'integral')", "Bubble(S, 4)", + "Bubble(S, 4, 'integral')", + "FacetBubble(T, 2)", + "FacetBubble(T, 2, 'integral')", + "FacetBubble(S, 3)", + "FacetBubble(S, 3, 'integral')", "RestrictedElement(Lagrange(I, 2), restriction_domain='facet')", "RestrictedElement(Lagrange(T, 2), restriction_domain='vertex')", "RestrictedElement(Lagrange(T, 3), restriction_domain='facet')", @@ -393,7 +400,6 @@ def __init__(self, a, b): xfail_impl("TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2))"), xfail_impl("Hdiv(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"), xfail_impl("Hcurl(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"), - xfail_impl("HDivTrace(T, 1)"), xfail_impl("EnrichedElement(" "Hdiv(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " "Hdiv(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" @@ -402,15 +408,16 @@ def __init__(self, a, b): "Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " "Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" ")"), + "HDivTrace(I, 0)", + "HDivTrace(T, 0)", + "HDivTrace(T, 1)", + "HDivTrace(T, 2)", + "HDivTrace(T, 3)", + "HDivTrace(S, 0)", + "HDivTrace(S, 1)", + "HDivTrace(S, 2)", + "HDivTrace(S, 3)", # Following elements are checked using tabulate - xfail_impl("HDivTrace(T, 0)"), - xfail_impl("HDivTrace(T, 1)"), - xfail_impl("HDivTrace(T, 2)"), - xfail_impl("HDivTrace(T, 3)"), - xfail_impl("HDivTrace(S, 0)"), - xfail_impl("HDivTrace(S, 1)"), - xfail_impl("HDivTrace(S, 2)"), - xfail_impl("HDivTrace(S, 3)"), xfail_impl("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))"), xfail_impl("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))"), xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))"), @@ -586,7 +593,6 @@ def test_facet_nodality_tabulate(element): facet_dim = sorted(entity_dofs.keys())[-2] facet_dofs = entity_dofs[facet_dim] dofs = element.dual_basis() - vertices = element.ref_el.vertices for (facet, indices) in facet_dofs.items(): for i in indices: @@ -595,16 +601,13 @@ def test_facet_nodality_tabulate(element): (coords, weights), = node.get_point_dict().items() assert weights == [(1.0, ())] - # Map dof coordinates to reference element due to - # HdivTrace interface peculiarity - ref_coords, = map_to_reference_facet((coords,), vertices, facet) - nodes_coords.append((facet, ref_coords)) + nodes_coords.append(coords) # Check nodality - for j, (facet, x) in enumerate(nodes_coords): - basis, = element.tabulate(0, (x,), entity=(facet_dim, facet)).values() - for i in range(len(basis)): - assert np.isclose(basis[i], 1.0 if i == j else 0.0) + expected = np.eye(len(dofs)) + for j, x in enumerate(nodes_coords): + basis, = element.tabulate(0, (x,)).values() + assert np.allclose(basis.T, expected[j]) @pytest.mark.parametrize('element', [ diff --git a/test/FIAT/unit/test_hdivtrace.py b/test/FIAT/unit/test_hdivtrace.py index ab66f0369..9c5be7505 100644 --- a/test/FIAT/unit/test_hdivtrace.py +++ b/test/FIAT/unit/test_hdivtrace.py @@ -25,7 +25,8 @@ @pytest.mark.parametrize("dim", (2, 3)) @pytest.mark.parametrize("degree", range(7)) -def test_basis_values(dim, degree): +@pytest.mark.parametrize("variant", ("spectral", "integral")) +def test_basis_values(dim, degree, variant): """Ensure that integrating simple monomials produces the expected results for each facet entity of the reference triangle and tetrahedron. @@ -35,12 +36,13 @@ def test_basis_values(dim, degree): (2) The entity pair (dim, id) is provided, and the trace element tabulates accordingly using the new tabulate API. """ - from FIAT import ufc_simplex, HDivTrace, make_quadrature + from FIAT import ufc_simplex, HDivTrace, make_quadrature, DiscontinuousLagrange ref_el = ufc_simplex(dim) - quadrule = make_quadrature(ufc_simplex(dim - 1), degree + 1) - fiat_element = HDivTrace(ref_el, degree) - facet_element = fiat_element.dg_elements[dim - 1] + facet = ref_el.get_facet_element() + quadrule = make_quadrature(facet, degree + 1) + fiat_element = HDivTrace(ref_el, degree, variant=variant) + facet_element = DiscontinuousLagrange(facet, degree, variant=variant) nf = facet_element.space_dimension() for facet_id in range(dim + 1): @@ -69,20 +71,22 @@ def test_basis_values(dim, degree): @pytest.mark.parametrize("degree", range(4)) -def test_quad_trace(degree): +@pytest.mark.parametrize("variant", ("spectral", "integral")) +def test_quad_trace(degree, variant): """Test the trace element defined on a quadrilateral cell""" - from FIAT import ufc_simplex, HDivTrace, make_quadrature - from FIAT.reference_element import TensorProductCell + from FIAT import ufc_cell, HDivTrace, make_quadrature, DiscontinuousLagrange + + ref_el = ufc_cell("quadrilateral") + fiat_element = HDivTrace(ref_el, degree, variant=variant) - tpc = TensorProductCell(ufc_simplex(1), ufc_simplex(1)) - fiat_element = HDivTrace(tpc, (degree, degree)) - facet_elements = fiat_element.dg_elements - quadrule = make_quadrature(ufc_simplex(1), degree + 1) + facet = ref_el.get_facet_element() + facet_element = DiscontinuousLagrange(facet, degree, variant=variant) + quadrule = make_quadrature(facet, degree + 1) - for i, entity in enumerate([((0, 1), 0), ((0, 1), 1), - ((1, 0), 0), ((1, 0), 1)]): - entity_dim, _ = entity - element = facet_elements[entity_dim] + entity_dim = facet.get_dimension() + for i, entity_id in enumerate(ref_el.topology[entity_dim]): + entity = (entity_dim, entity_id) + element = facet_element nf = element.space_dimension() tab = fiat_element.tabulate(0, quadrule.pts,