From c073ac56ef9ce1a8c6b87543d5fa79b83331c1dc Mon Sep 17 00:00:00 2001 From: Jack Betteridge <43041811+JDBetteridge@users.noreply.github.com> Date: Wed, 4 Sep 2024 16:21:27 +0100 Subject: [PATCH 01/13] UPDATE Python (#81) * Setuptools is now required to provide pkg_resources module * Update python versions used for CI --- .github/workflows/pythonapp.yml | 18 +++++++++--------- setup.py | 24 ++++++++++++++---------- 2 files changed, 23 insertions(+), 19 deletions(-) diff --git a/.github/workflows/pythonapp.yml b/.github/workflows/pythonapp.yml index 39dd4bcc9..ce154290f 100644 --- a/.github/workflows/pythonapp.yml +++ b/.github/workflows/pythonapp.yml @@ -12,30 +12,30 @@ jobs: strategy: matrix: os: [ubuntu-latest, macos-latest] - python-version: [3.8, 3.9] + python-version: ['3.9', '3.10', '3.11', '3.12'] steps: - - uses: actions/checkout@v2 + - uses: actions/checkout@v4 - name: Set up Python - uses: actions/setup-python@v2 + uses: actions/setup-python@v5 with: python-version: ${{ matrix.python-version }} - name: Lint with flake8 run: | - pip install flake8 - flake8 --statistics . + python -m pip install flake8 + python -m flake8 --statistics . - name: Check documentation style run: | - pip install pydocstyle + python -m pip install pydocstyle python -m pydocstyle . - name: Install FIAT run: pip install . - name: Test with pytest run: | - pip install coveralls pytest pytest-cov pytest-xdist - DATA_REPO_GIT="" pytest --cov=FIAT/ test/ + python -m pip install coveralls pytest pytest-cov pytest-xdist + DATA_REPO_GIT="" python -m pytest --cov=FIAT/ test/ - name: Coveralls - if: ${{ github.repository == 'FEniCS/fiat' && github.head_ref == '' && matrix.os == 'ubuntu-latest' && matrix.python-version == '3.8' }} + if: ${{ github.repository == 'FEniCS/fiat' && github.head_ref == '' && matrix.os == 'ubuntu-latest' && matrix.python-version == '3.11' }} env: COVERALLS_REPO_TOKEN: ${{ secrets.COVERALLS_REPO_TOKEN }} run: coveralls diff --git a/setup.py b/setup.py index 20ba60391..d55401a90 100755 --- a/setup.py +++ b/setup.py @@ -18,13 +18,17 @@ if 'dev' not in version: tarball = url + "downloads/fenics-fiat-%s.tar.gz" % version -setup(name="fenics-fiat", - description="FInite element Automatic Tabulator", - version=version, - author="Robert C. Kirby et al.", - author_email="fenics-dev@googlegroups.com", - url=url, - download_url=tarball, - license="LGPL v3 or later", - packages=["FIAT"], - install_requires=["numpy", "recursivenodes", "scipy", "sympy"]) +setup( + name="fenics-fiat", + description="FInite element Automatic Tabulator", + version=version, + author="Robert C. Kirby et al.", + author_email="fenics-dev@googlegroups.com", + url=url, + download_url=tarball, + license="LGPL v3 or later", + packages=["FIAT"], + install_requires=[ + "setuptools", "numpy", "recursivenodes", "scipy", "sympy" + ] +) From 78fc962a3b9d475cfd2ceb3000cc10028ed3c5bb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 11 Sep 2024 07:50:32 +0100 Subject: [PATCH 02/13] Macro: bin points to nearest cell (#82) * Macro: bin points to nearest cell * Compute l1 distances with the correct scale * Test compute_l1_distance * Generalize distance_to_point_l1 * clean up reference_element.py * optionally rescale the L1 distance * generalize compute_barycentric_coordinates * vectorize entity_transform --- FIAT/Sminus.py | 2 +- FIAT/SminusCurl.py | 2 +- FIAT/SminusDiv.py | 2 +- FIAT/bernstein.py | 2 +- FIAT/dual_set.py | 2 +- FIAT/expansions.py | 76 ++--- FIAT/finite_element.py | 4 +- FIAT/reference_element.py | 343 +++++++++++------------ FIAT/serendipity.py | 3 +- test/unit/test_gauss_legendre.py | 2 +- test/unit/test_gauss_lobatto_legendre.py | 2 +- test/unit/test_hdivtrace.py | 2 +- test/unit/test_macro.py | 33 ++- test/unit/test_reference_element.py | 3 +- 14 files changed, 237 insertions(+), 241 deletions(-) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index a7e7735ac..6d9d4eaa3 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -150,7 +150,7 @@ def tabulate(self, order, points, entity=None): entity_dim, entity_id = entity transform = self.ref_el.get_entity_transform(entity_dim, entity_id) - points = list(map(transform, points)) + points = transform(points) phivals = {} diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index e6fd93d06..ce42a194f 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -131,7 +131,7 @@ def tabulate(self, order, points, entity=None): entity_dim, entity_id = entity transform = self.ref_el.get_entity_transform(entity_dim, entity_id) - points = list(map(transform, points)) + points = transform(points) phivals = {} diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 7ccbc8aae..70527a09b 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -121,7 +121,7 @@ def tabulate(self, order, points, entity=None): entity_dim, entity_id = entity transform = self.ref_el.get_entity_transform(entity_dim, entity_id) - points = list(map(transform, points)) + points = transform(points) phivals = {} diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index 7fe6cfdb4..81768d12d 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -80,7 +80,7 @@ def tabulate(self, order, points, entity=None): entity_dim, entity_id = entity entity_transform = ref_el.get_entity_transform(entity_dim, entity_id) - cell_points = list(map(entity_transform, points)) + cell_points = entity_transform(points) # Construct Cartesian to Barycentric coordinate mapping vs = numpy.asarray(ref_el.get_vertices()) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index bf038ac62..d6d0de44b 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -258,7 +258,7 @@ def lexsort_nodes(ref_el, nodes, entity=None, offset=0): for node in nodes: pt, = node.get_point_dict() pts.append(pt) - bary = ref_el.compute_barycentric_coordinates(pts, entity=entity) + bary = ref_el.compute_barycentric_coordinates(pts) order = list(offset + numpy.lexsort(bary.T)) else: order = list(range(offset, offset + len(nodes))) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index bf38f9795..8ff65f164 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -242,21 +242,6 @@ def xi_tetrahedron(eta): return xi1, xi2, xi3 -def apply_mapping(A, b, pts, transpose=False): - """Apply an affine mapping to a column-stacked array of points.""" - if len(pts) == 0: - return pts - if transpose: - Ax = numpy.dot(pts, A.T) - for _ in Ax.shape[:-1]: - b = b[None, ...] - else: - Ax = numpy.dot(A, pts) - for _ in Ax.shape[1:]: - b = b[..., None] - return Ax + b - - class ExpansionSet(object): def __new__(cls, *args, **kwargs): """Returns an ExpansionSet instance appropriate for the given @@ -322,7 +307,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): from FIAT.polynomial_set import mis lorder = min(order, self.recurrence_order) A, b = self.affine_mappings[cell] - ref_pts = apply_mapping(A, b, numpy.transpose(pts)) + ref_pts = numpy.add(numpy.dot(pts, A.T), b).T Jinv = A if direction is None else numpy.dot(A, direction)[:, None] sd = self.ref_el.get_spatial_dimension() @@ -411,7 +396,7 @@ def tabulate_normal_jumps(self, n, ref_pts, facet, order=0): """ sd = self.ref_el.get_spatial_dimension() transform = self.ref_el.get_entity_transform(sd-1, facet) - pts = numpy.array(list(map(transform, ref_pts))) + pts = transform(ref_pts) cell_point_map = compute_cell_point_map(self.ref_el, pts, unique=False) cell_node_map = self.get_cell_node_map(n) @@ -557,7 +542,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): A, b = self.affine_mappings[cell] Jinv = A[0, 0] if direction is None else numpy.dot(A, direction) - xs = apply_mapping(A, b, numpy.transpose(pts)).T + xs = numpy.add(numpy.dot(pts, A.T), b) results = {} scale = self.get_scale(cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1) for k in range(order+1): @@ -656,39 +641,15 @@ def polynomial_cell_node_map(ref_el, n, continuity=None): return cell_node_map -def compute_l1_distance(ref_el, points, entity=None): - """Computes the l1 distances from a simplex entity to a set of points. - - :arg ref_el: a SimplicialComplex. - :arg points: an iterable of points. - :arg entity: a tuple of the entity dimension and id. - :returns: a numpy array with the l1 distance from the entity to each point. - """ - if entity is None: - entity = (ref_el.get_spatial_dimension(), 0) - dim, entity_id = entity - ref_verts = numpy.zeros((dim + 1, dim)) - ref_verts[range(1, dim + 1), range(dim)] = 1.0 - - top = ref_el.get_topology() - verts = ref_el.get_vertices_of_subcomplex(top[dim][entity_id]) - A, b = reference_element.make_affine_mapping(verts, ref_verts) - A = numpy.vstack((A, -numpy.sum(A, axis=0))) - b = numpy.hstack((b, 1-numpy.sum(b, axis=0))) - bary = apply_mapping(A, b, points, transpose=True) - - dist = 0.5 * abs(numpy.sum(abs(bary) - bary, axis=-1)) - return dist - - def compute_cell_point_map(ref_el, pts, unique=True, tol=1E-12): """Maps cells on a simplicial complex to points. + Points outside the complex are binned to the nearest cell. :arg ref_el: a SimplicialComplex. :arg pts: an iterable of physical points on the complex. :kwarg unique: Are we assigning a unique cell to points on facets? :kwarg tol: the absolute tolerance. - :returns: a dict mapping cell id to points located on that cell. + :returns: a dict mapping cell id to the point ids nearest to that cell. """ top = ref_el.get_topology() sd = ref_el.get_spatial_dimension() @@ -699,21 +660,25 @@ def compute_cell_point_map(ref_el, pts, unique=True, tol=1E-12): if pts.dtype == object: 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) + tol = best + tol + cell_point_map = {} for cell in sorted(top[sd]): # Bin points based on l1 distance - pts_on_cell = compute_l1_distance(ref_el, pts, entity=(sd, cell)) < tol - if len(pts_on_cell.shape) == 0: + pts_near_cell = ref_el.distance_to_point_l1(pts, entity=(sd, cell), rescale=True) < tol + if len(pts_near_cell.shape) == 0: # singleton case - if pts_on_cell: + if pts_near_cell: cell_point_map[cell] = Ellipsis if unique: break else: if unique: for other in cell_point_map.values(): - pts_on_cell[other] = False - ipts = numpy.where(pts_on_cell)[0] + pts_near_cell[other] = False + ipts = numpy.where(pts_near_cell)[0] if len(ipts) > 0: cell_point_map[cell] = ipts return cell_point_map @@ -734,14 +699,19 @@ def compute_partition_of_unity(ref_el, pt, unique=True, tol=1E-12): # assert singleton point pt = pt.reshape((sd,)) - # Compute characteristic function of each cell + # The distance to the nearest cell is equal to the distance to the parent cell + best = ref_el.get_parent().distance_to_point_l1(pt, rescale=True) + tol = best + tol + + # Compute characteristic function of each subcell otherwise = [] masks = [] for cell in sorted(top[sd]): - inside = compute_l1_distance(ref_el, pt, entity=(sd, cell)) < tol - masks.append(Piecewise(*otherwise, (1.0, inside), (0.0, True))) + # Bin points based on l1 distance + pt_near_cell = ref_el.distance_to_point_l1(pt, entity=(sd, cell), rescale=True) < tol + masks.append(Piecewise(*otherwise, (1.0, pt_near_cell), (0.0, True))) if unique: - otherwise.append((0.0, inside)) + otherwise.append((0.0, pt_near_cell)) # If the point is on a facet, divide the characteristic function by the facet multiplicity if not unique: mult = sum(masks) diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index 40b2bcb1c..20d05f142 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -191,7 +191,7 @@ def tabulate(self, order, points, entity=None): entity_dim, entity_id = entity transform = self.ref_el.get_entity_transform(entity_dim, entity_id) - return self.poly_set.tabulate(list(map(transform, points)), order) + return self.poly_set.tabulate(transform(points), order) def value_shape(self): "Return the value shape of the finite element functions." @@ -243,7 +243,7 @@ def entity_support_dofs(elem, entity_dim): result = {} for f in elem.entity_dofs()[entity_dim].keys(): entity_transform = ref_el.get_entity_transform(entity_dim, f) - points = list(map(entity_transform, quad.get_points())) + points = entity_transform(quad.get_points()) # Integrate the square of the basis functions on the facet. vals = elem.tabulate(0, points)[(0,) * dim] diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 379b2958a..1e750bbed 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -377,9 +377,8 @@ def compute_tangents(self, dim, i): of dimension dim. Returns a (possibly empty) list. These tangents are *NOT* normalized to have unit length.""" t = self.get_topology() - vs = list(map(numpy.array, self.get_vertices_of_subcomplex(t[dim][i]))) - ts = [v - vs[0] for v in vs[1:]] - return ts + vs = numpy.array(self.get_vertices_of_subcomplex(t[dim][i])) + return vs[1:] - vs[:1] def compute_normalized_tangents(self, dim, i): """Computes tangents in any dimension based on differences @@ -387,19 +386,21 @@ def compute_normalized_tangents(self, dim, i): of dimension dim. Returns a (possibly empty) list. These tangents are normalized to have unit length.""" ts = self.compute_tangents(dim, i) - return [t / numpy.linalg.norm(t) for t in ts] + ts /= numpy.linalg.norm(ts, axis=0)[None, :] + return ts def compute_edge_tangent(self, edge_i): """Computes the nonnormalized tangent to a 1-dimensional facet. returns a single vector.""" t = self.get_topology() - (v0, v1) = self.get_vertices_of_subcomplex(t[1][edge_i]) - return numpy.array(v1) - numpy.array(v0) + vs = numpy.asarray(self.get_vertices_of_subcomplex(t[1][edge_i])) + return vs[1] - vs[0] def compute_normalized_edge_tangent(self, edge_i): """Computes the unit tangent vector to a 1-dimensional facet""" v = self.compute_edge_tangent(edge_i) - return v / numpy.linalg.norm(v) + v /= numpy.linalg.norm(v) + return v def compute_face_tangents(self, face_i): """Computes the two tangents to a face. Only implemented @@ -407,9 +408,8 @@ def compute_face_tangents(self, face_i): if self.get_spatial_dimension() != 3: raise Exception("can't get face tangents yet") t = self.get_topology() - (v0, v1, v2) = list(map(numpy.array, - self.get_vertices_of_subcomplex(t[2][face_i]))) - return (v1 - v0, v2 - v0) + vs = numpy.asarray(self.get_vertices_of_subcomplex(t[2][face_i])) + return vs[1:] - vs[:1] def compute_face_edge_tangents(self, dim, entity_id): """Computes all the edge tangents of any k-face with k>=1. @@ -417,13 +417,14 @@ def compute_face_edge_tangents(self, dim, entity_id): This agrees with `compute_edge_tangent` when dim=1. """ vert_ids = self.get_topology()[dim][entity_id] - vert_coords = [numpy.array(x) - for x in self.get_vertices_of_subcomplex(vert_ids)] - edge_ts = [] + vert_coords = numpy.asarray(self.get_vertices_of_subcomplex(vert_ids)) + v0 = [] + v1 = [] for source in range(dim): for dest in range(source + 1, dim + 1): - edge_ts.append(vert_coords[dest] - vert_coords[source]) - return edge_ts + v0.append(source) + v1.append(dest) + return vert_coords[v1] - vert_coords[v0] def make_points(self, dim, entity_id, order, variant=None): """Constructs a lattice of points on the entity_id:th @@ -476,156 +477,87 @@ def get_entity_transform(self, dim, entity): if dim == 0: # Special case vertices. i, = topology[dim][entity] - vertex = self.get_vertices()[i] - return lambda point: vertex + offset = numpy.asarray(self.get_vertices()[i]) + C = numpy.zeros((dim, ) + offset.shape) elif dim == celldim and len(self.topology[celldim]) == 1: assert entity == 0 - return lambda point: point - - try: + return lambda x: x + else: subcell = self.construct_subelement(dim) - except NotImplementedError: - # Special case for 1D elements. - x_c, = self.get_vertices_of_subcomplex(topology[0][entity]) - return lambda x: x_c - - subdim = subcell.get_spatial_dimension() + subdim = subcell.get_spatial_dimension() + assert subdim == celldim - codim - assert subdim == celldim - codim + # Entity vertices in entity space. + v_e = numpy.asarray(subcell.get_vertices()) + A = v_e[1:] - v_e[:1] - # Entity vertices in entity space. - v_e = numpy.asarray(subcell.get_vertices()) - A = v_e[1:] - v_e[:1] + # Entity vertices in cell space. + v_c = numpy.asarray(self.get_vertices_of_subcomplex(topology[dim][entity])) + B = v_c[1:] - v_c[:1] - # Entity vertices in cell space. - v_c = numpy.asarray(self.get_vertices_of_subcomplex(topology[dim][entity])) - B = v_c[1:] - v_c[:1] + C = numpy.linalg.solve(A, B) - C = numpy.linalg.solve(A, B).T + offset = v_c[0] - numpy.dot(v_e[0], C) - offset = v_c[0] - C.dot(v_e[0]) + def transform(point): + out = numpy.dot(point, C) + return numpy.add(out, offset, out=out) - return lambda x: offset + C.dot(x) + return transform def get_dimension(self): """Returns the subelement dimension of the cell. Same as the spatial dimension.""" return self.get_spatial_dimension() - def compute_barycentric_coordinates(self, points, entity=None): + def compute_barycentric_coordinates(self, points, entity=None, rescale=False): """Returns the barycentric coordinates of a list of points on an entity.""" + if len(points) == 0: + return points if entity is None: entity = (self.get_spatial_dimension(), 0) entity_dim, entity_id = entity top = self.get_topology() verts = self.get_vertices_of_subcomplex(top[entity_dim][entity_id]) - A = numpy.transpose(verts) - B = numpy.transpose(points) - B = B - A[:, :1] - A = A[:, 1:] - A[:, :1] - if A.shape[0] != A.shape[1]: - # Form normal equations - B = numpy.dot(A.T, B) - A = numpy.dot(A.T, A) - return numpy.linalg.solve(A, B).T - - -class Simplex(SimplicialComplex): - r"""Abstract class for a reference simplex. - - Orientation of a physical cell is computed systematically - by comparing the canonical orderings of its facets and - the facets in the FIAT reference cell. - - As an example, we compute the orientation of a - triangular cell: - - + + - | \ | \ - 1 0 47 42 - | \ | \ - +--2---+ +--43--+ - FIAT canonical Mapped example physical cell - - Suppose that the facets of the physical cell - are canonically ordered as: - - C = [43, 42, 47] - - FIAT facet to Physical facet map is given by: - - M = [42, 47, 43] - - Then the orientation of the cell is computed as: - - C.index(M[0]) = 1; C.remove(M[0]) - C.index(M[1]) = 1; C.remove(M[1]) - C.index(M[2]) = 0; C.remove(M[2]) - - o = (1 * 2!) + (1 * 1!) + (0 * 0!) = 3 - """ - def is_simplex(self): - return True - - def symmetry_group_size(self, dim): - return factorial(dim + 1) - - def cell_orientation_reflection_map(self): - """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" - return make_cell_orientation_reflection_map_simplex(self.get_dimension()) - - -# Backwards compatible name -ReferenceElement = Simplex - - -class UFCSimplex(Simplex): - - def get_facet_element(self): - dimension = self.get_spatial_dimension() - return self.construct_subelement(dimension - 1) - - def construct_subelement(self, dimension): - """Constructs the reference element of a cell subentity - specified by subelement dimension. - - :arg dimension: subentity dimension (integer) - """ - return ufc_simplex(dimension) - - def contains_point(self, point, epsilon=0): - """Checks if reference cell contains given point - (with numerical tolerance as given by the L1 distance (aka 'manhatten', - 'taxicab' or rectilinear distance) to the cell. - - Parameters - ---------- - point : numpy.ndarray, list or symbolic expression - The coordinates of the point. - epsilon : float - The tolerance for the check. - - Returns - ------- - bool : True if the point is inside the cell, False otherwise. - - """ - return self.distance_to_point_l1(point) <= epsilon - - def distance_to_point_l1(self, point): + if entity_dim == self.get_spatial_dimension(): + ref_verts = numpy.eye(entity_dim + 1) + A, b = make_affine_mapping(verts, ref_verts) + else: + v = numpy.transpose(verts) + v = v[:, 1:] - v[:, :1] + A = numpy.linalg.solve(numpy.dot(v.T, v), v.T) + b = -numpy.dot(A, verts[0]) + A = numpy.vstack((-numpy.sum(A, axis=0), A)) + b = numpy.hstack((1 - numpy.sum(b, axis=0), b)) + + if rescale: + # rescale barycentric coordinates by the height wrt. to 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) + + def distance_to_point_l1(self, points, entity=None, rescale=False): # noqa: D301 """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear - distance) to a point with 0.0 if the point is inside the cell. + distance) from an entity to a point with 0.0 if the point is inside the entity. Parameters ---------- - point : numpy.ndarray or list - The coordinates of the point. + points : numpy.ndarray or list + The coordinates of the points. + entity : tuple or None + A tuple of entity dimension and entity id. + rescale : bool + If true, the L1 distance is measured with respect to rescaled + barycentric coordinates, such that the L1 and L2 distances agree + for points opposite to a single facet. Returns ------- - float + numpy.float64 or numpy.ndarray The L1 distance, also known as taxicab, manhatten or rectilinear distance, of the cell to the point. If 0.0 the point is inside the cell. @@ -728,35 +660,104 @@ def distance_to_point_l1(self, point): `beta = X[0] = x`, `gamma = X[1] = y` and `delta = X[2] = z`. - The rules are the same as for the triangle but with one extra + The rules are the same as for the tetrahedron but with one extra barycentric coordinate. Our approximate distance, the absolute sum of the negative barycentric coordinates, is at worse around 4 times the actual distance to the tetrahedron. """ - # bary = [alpha, beta, gamma, delta, ...] - see docstring - bary = [1.0 - sum(point)] + list(point) - # We avoid branching so that code can be generated from this (e.g. with - # sympy). bary-abs(bary) gets rid of the positive barycentric - # coordinates and doubles the negative distances. Summing, halving and - # taking the negative of these gives the L1 distance. So for example - # point = [-1, -1] - # bary = [3, -1, -1], - # bary-abs(bary) = [0, -2, -2], - # sum(bary-abs(bary)) = -4. - # - 0.5 * sum(bary-abs(bary)) = 2.0 - # which is the correct L1 distance from the cell to the point. - l1_dist = - 0.5 * sum(b - abs(b) for b in bary) - # Take abs at the end to avoid negative zero. - return abs(l1_dist) + # sum the negative part of each barycentric coordinate + bary = self.compute_barycentric_coordinates(points, entity=entity, rescale=rescale) + return 0.5 * abs(numpy.sum(abs(bary) - bary, axis=-1)) + def contains_point(self, point, epsilon=0.0, entity=None): + """Checks if reference cell contains given point + (with numerical tolerance as given by the L1 distance (aka 'manhatten', + 'taxicab' or rectilinear distance) to the cell. -class DefaultSimplex(Simplex): + Parameters + ---------- + point : numpy.ndarray, list or symbolic expression + The coordinates of the point. + epsilon : float + The tolerance for the check. + entity : tuple or None + A tuple of entity dimension and entity id. + + Returns + ------- + bool : True if the point is inside the cell, False otherwise. + + """ + return self.distance_to_point_l1(point, entity=entity) <= epsilon + + +class Simplex(SimplicialComplex): + r"""Abstract class for a reference simplex. + + Orientation of a physical cell is computed systematically + by comparing the canonical orderings of its facets and + the facets in the FIAT reference cell. + + As an example, we compute the orientation of a + triangular cell: + + + + + | \ | \ + 1 0 47 42 + | \ | \ + +--2---+ +--43--+ + FIAT canonical Mapped example physical cell + + Suppose that the facets of the physical cell + are canonically ordered as: + + C = [43, 42, 47] + + FIAT facet to Physical facet map is given by: + + M = [42, 47, 43] + + Then the orientation of the cell is computed as: + + C.index(M[0]) = 1; C.remove(M[0]) + C.index(M[1]) = 1; C.remove(M[1]) + C.index(M[2]) = 0; C.remove(M[2]) + + o = (1 * 2!) + (1 * 1!) + (0 * 0!) = 3 + """ + def is_simplex(self): + return True + + def symmetry_group_size(self, dim): + return factorial(dim + 1) + + def cell_orientation_reflection_map(self): + """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" + return make_cell_orientation_reflection_map_simplex(self.get_dimension()) def get_facet_element(self): dimension = self.get_spatial_dimension() return self.construct_subelement(dimension - 1) + +# Backwards compatible name +ReferenceElement = Simplex + + +class UFCSimplex(Simplex): + + def construct_subelement(self, dimension): + """Constructs the reference element of a cell subentity + specified by subelement dimension. + + :arg dimension: subentity dimension (integer) + """ + return ufc_simplex(dimension) + + +class DefaultSimplex(Simplex): + def construct_subelement(self, dimension): """Constructs the reference element of a cell subentity specified by subelement dimension. @@ -768,10 +769,6 @@ def construct_subelement(self, dimension): class SymmetricSimplex(Simplex): - def get_facet_element(self): - dimension = self.get_spatial_dimension() - return self.construct_subelement(dimension - 1) - def construct_subelement(self, dimension): """Constructs the reference element of a cell subentity specified by subelement dimension. @@ -1042,8 +1039,10 @@ def get_entity_transform(self, dim, entity_i): slices = TensorProductCell._split_slices(dim) def transform(point): - return list(chain(*[t(point[s]) - for t, s in zip(sct, slices)])) + point = numpy.asarray(point) + return numpy.concatenate(tuple(t(point[..., s]) + for t, s in zip(sct, slices)), axis=-1) + return transform def volume(self): @@ -1065,7 +1064,7 @@ def compute_reference_normal(self, facet_dim, facet_i): n.extend([0] * c.get_spatial_dimension()) return numpy.asarray(n) - def contains_point(self, point, epsilon=0): + def contains_point(self, point, epsilon=0.0): """Checks if reference cell contains given point (with numerical tolerance as given by the L1 distance (aka 'manhatten', 'taxicab' or rectilinear distance) to the cell. @@ -1082,27 +1081,25 @@ def contains_point(self, point, epsilon=0): bool : True if the point is inside the cell, False otherwise. """ - subcell_dimensions = [c.get_spatial_dimension() for c in self.cells] + subcell_dimensions = self.get_dimension() assert len(point) == sum(subcell_dimensions) point_slices = TensorProductCell._split_slices(subcell_dimensions) - subcell_points = [point[s] for s in point_slices] return reduce(operator.and_, - (c.contains_point(p, epsilon=epsilon) - for c, p in zip(self.cells, subcell_points)), + (c.contains_point(point[s], epsilon=epsilon) + for c, s in zip(self.cells, point_slices)), True) - def distance_to_point_l1(self, point): + def distance_to_point_l1(self, point, rescale=False): """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear distance) to a point with 0.0 if the point is inside the cell. For more information see the docstring for the UFCSimplex method.""" - subcell_dimensions = [c.get_spatial_dimension() for c in self.cells] + subcell_dimensions = self.get_dimension() assert len(point) == sum(subcell_dimensions) point_slices = TensorProductCell._split_slices(subcell_dimensions) - subcell_points = [point[s] for s in point_slices] - subcell_distances = [c.distance_to_point_l1(p) - for c, p in zip(self.cells, subcell_points)] - return sum(subcell_distances) + point = numpy.asarray(point) + return sum(c.distance_to_point_l1(point[..., s], rescale=rescale) + for c, s in zip(self.cells, point_slices)) def symmetry_group_size(self, dim): return tuple(c.symmetry_group_size(d) for d, c in zip(dim, self.cells)) @@ -1248,12 +1245,12 @@ def contains_point(self, point, epsilon=0): """ return self.product.contains_point(point, epsilon=epsilon) - def distance_to_point_l1(self, point): + def distance_to_point_l1(self, point, rescale=False): """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear distance) to a point with 0.0 if the point is inside the cell. For more information see the docstring for the UFCSimplex method.""" - return self.product.distance_to_point_l1(point) + return self.product.distance_to_point_l1(point, rescale=rescale) def symmetry_group_size(self, dim): return [1, 2, 8][dim] @@ -1353,12 +1350,12 @@ def contains_point(self, point, epsilon=0): """ return self.product.contains_point(point, epsilon=epsilon) - def distance_to_point_l1(self, point): + def distance_to_point_l1(self, point, rescale=False): """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear distance) to a point with 0.0 if the point is inside the cell. For more information see the docstring for the UFCSimplex method.""" - return self.product.distance_to_point_l1(point) + return self.product.distance_to_point_l1(point, rescale=rescale) def symmetry_group_size(self, dim): return [1, 2, 8, 48][dim] diff --git a/FIAT/serendipity.py b/FIAT/serendipity.py index 30314519f..4f783dd58 100644 --- a/FIAT/serendipity.py +++ b/FIAT/serendipity.py @@ -169,7 +169,7 @@ def tabulate(self, order, points, entity=None): entity_dim, entity_id = entity transform = self.ref_el.get_entity_transform(entity_dim, entity_id) - points = list(map(transform, points)) + points = transform(points) phivals = {} dim = self.flat_el.get_spatial_dimension() @@ -177,7 +177,6 @@ def tabulate(self, order, points, entity=None): raise NotImplementedError('no tabulate method for serendipity elements of dimension 1 or less.') if dim >= 4: raise NotImplementedError('tabulate does not support higher dimensions than 3.') - points = np.asarray(points) npoints, pointdim = points.shape for o in range(order + 1): alphas = mis(dim, o) diff --git a/test/unit/test_gauss_legendre.py b/test/unit/test_gauss_legendre.py index 676a4ea26..36fd53bc9 100644 --- a/test/unit/test_gauss_legendre.py +++ b/test/unit/test_gauss_legendre.py @@ -66,7 +66,7 @@ def test_edge_dofs(dim, degree): for entity in edge_dofs: if len(edge_dofs[entity]) > 0: transform = s.get_entity_transform(1, entity) - assert np.allclose(points[edge_dofs[entity]], np.array(list(map(transform, quadrature_points)))) + assert np.allclose(points[edge_dofs[entity]], transform(quadrature_points)) @pytest.mark.parametrize("dim, degree", [(1, 64), (2, 16), (3, 16)]) diff --git a/test/unit/test_gauss_lobatto_legendre.py b/test/unit/test_gauss_lobatto_legendre.py index 4c288d654..5a7d0e07d 100644 --- a/test/unit/test_gauss_lobatto_legendre.py +++ b/test/unit/test_gauss_lobatto_legendre.py @@ -66,7 +66,7 @@ def test_edge_dofs(dim, degree): for entity in edge_dofs: if len(edge_dofs[entity]) > 0: transform = s.get_entity_transform(1, entity) - assert np.allclose(points[edge_dofs[entity]], np.array(list(map(transform, quadrature_points)))) + assert np.allclose(points[edge_dofs[entity]], transform(quadrature_points)) @pytest.mark.parametrize("dim, degree", [(1, 64), (2, 16), (3, 16)]) diff --git a/test/unit/test_hdivtrace.py b/test/unit/test_hdivtrace.py index d85c47d62..ab66f0369 100644 --- a/test/unit/test_hdivtrace.py +++ b/test/unit/test_hdivtrace.py @@ -46,7 +46,7 @@ def test_basis_values(dim, degree): for facet_id in range(dim + 1): # Tabulate without an entity pair given --- need to map to cell coordinates cell_transform = ref_el.get_entity_transform(dim - 1, facet_id) - cell_points = np.array(list(map(cell_transform, quadrule.pts))) + cell_points = cell_transform(quadrule.get_points()) ctab = fiat_element.tabulate(0, cell_points)[(0,) * dim][nf*facet_id:nf*(facet_id + 1)] # Tabulate with entity pair provided diff --git a/test/unit/test_macro.py b/test/unit/test_macro.py index cc0dada0f..4a6646716 100644 --- a/test/unit/test_macro.py +++ b/test/unit/test_macro.py @@ -43,7 +43,7 @@ def test_split_make_points(split, cell, degree, variant): for entity in top[i]: pts_entity = split_cell.make_points(i, entity, degree, variant=variant) mapping = split_cell.get_entity_transform(i, entity) - mapped_pts = list(map(mapping, pts_ref)) + mapped_pts = mapping(pts_ref) assert numpy.allclose(mapped_pts, pts_entity) @@ -343,6 +343,37 @@ def test_Ck_basis(cell, order, degree, variant): assert numpy.allclose(local_phis, phis[:, ipts]) +def test_distance_to_point_l1(cell): + A = AlfeldSplit(cell) + dim = A.get_spatial_dimension() + top = A.get_topology() + p0, = cell.make_points(dim, 0, dim+1) + + # construct one point in front of each facet + pts = [] + expected = [] + parent_top = cell.get_topology() + for i in parent_top[dim-1]: + Fi, = numpy.asarray(cell.make_points(dim-1, i, dim)) + n = cell.compute_normal(i) + n *= numpy.dot(n, Fi - p0) + n /= numpy.linalg.norm(n) + d = 0.222 + i/10 + pts.append(Fi + d * n) + expected.append(d) + + # the computed L1 distance agrees with the L2 distance for points in front of facets + parent_distance = cell.distance_to_point_l1(pts, rescale=True) + assert numpy.allclose(parent_distance, expected) + + # assert that the subcell measures the same distance as the parent + for i in top[dim]: + subcell_distance = A.distance_to_point_l1(pts, entity=(dim, i), rescale=True) + assert numpy.isclose(subcell_distance[i], expected[i]) + assert all(subcell_distance[:i] > expected[:i]) + assert all(subcell_distance[i+1:] > expected[i+1:]) + + @pytest.mark.parametrize("element", (DiscontinuousLagrange, Lagrange)) def test_macro_sympy(cell, element): import sympy diff --git a/test/unit/test_reference_element.py b/test/unit/test_reference_element.py index 2945508b2..f551f4d69 100644 --- a/test/unit/test_reference_element.py +++ b/test/unit/test_reference_element.py @@ -18,7 +18,6 @@ import pytest import numpy as np import sys -from math import isclose from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron from FIAT.reference_element import Point, TensorProductCell, UFCQuadrilateral, UFCHexahedron @@ -394,7 +393,7 @@ def test_contains_point(cell, point, epsilon, expected): (quadrilateral_x_interval, [0.0, 0.0, -1e-12], 1e-12), (quadrilateral_x_interval, [0.0, 0.0, 1+1e-12], 1e-12)]) def test_distance_to_point_l1(cell, point, expected): - assert isclose(cell.distance_to_point_l1(point), expected, rel_tol=1e-3) + assert np.isclose(cell.distance_to_point_l1(point), expected, rtol=1e-3) if __name__ == '__main__': From 13131da108396ae6ce6ef426c99547041682d427 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 8 Oct 2024 09:25:33 +0100 Subject: [PATCH 03/13] Add Stokes macroelements (#84) * Alfeld-Sorokina macroelement * Add the Christiansen-Hu macroelement (on barycentric Worsey-Farin splits) * PowellSabinSplit has AlfeldSplit as parent_complex * Add ArnoldQin / ReducedArnoldQin * add Bernardi-Raugel * clean up functionals.py * JM: fix scaling * AlfeldSplit and WorseyFarinSplit subclass PowellSabinSplit * Fix super * Powell-Sabin: Preserve the numbering of the unsplit entities * Add "worsey-farin" Lagrange variant * Add Guzman-Neilan * RestrictedElement: preserve ordering --- FIAT/P0.py | 4 +- FIAT/Sminus.py | 4 +- FIAT/SminusCurl.py | 12 +- FIAT/SminusDiv.py | 14 +- FIAT/__init__.py | 15 +- FIAT/alfeld_sorokina.py | 106 ++++++++++++++ FIAT/argyris.py | 4 +- FIAT/arnold_qin.py | 71 +++++++++ FIAT/arnold_winther.py | 17 +-- FIAT/barycentric_interpolation.py | 5 +- FIAT/bell.py | 6 +- FIAT/bernardi_raugel.py | 102 +++++++++++++ FIAT/bernstein.py | 4 +- FIAT/brezzi_douglas_fortin_marini.py | 6 +- FIAT/brezzi_douglas_marini.py | 5 +- FIAT/brezzi_douglas_marini_cube.py | 12 +- FIAT/bubble.py | 6 +- FIAT/check_format_variant.py | 18 ++- FIAT/christiansen_hu.py | 78 ++++++++++ FIAT/crouzeix_raviart.py | 4 +- FIAT/discontinuous_lagrange.py | 8 +- FIAT/discontinuous_pc.py | 22 +-- FIAT/discontinuous_raviart_thomas.py | 5 +- FIAT/discontinuous_taylor.py | 4 +- FIAT/enriched.py | 2 +- FIAT/expansions.py | 12 +- FIAT/fdm_element.py | 4 +- FIAT/functional.py | 212 ++++++++++++--------------- FIAT/gauss_legendre.py | 2 +- FIAT/gauss_lobatto_legendre.py | 2 +- FIAT/gauss_radau.py | 4 +- FIAT/guzman_neilan.py | 147 +++++++++++++++++++ FIAT/hct.py | 17 +-- FIAT/hdiv_trace.py | 8 +- FIAT/hellan_herrmann_johnson.py | 5 +- FIAT/hermite.py | 4 +- FIAT/hierarchical.py | 10 +- FIAT/johnson_mercier.py | 17 +-- FIAT/lagrange.py | 4 +- FIAT/macro.py | 160 ++++++++++---------- FIAT/nedelec.py | 5 +- FIAT/nedelec_second_kind.py | 4 +- FIAT/powell_sabin.py | 10 +- FIAT/quadrature.py | 12 +- FIAT/quadrature_element.py | 2 +- FIAT/raviart_thomas.py | 7 +- FIAT/reference_element.py | 82 ++++++++--- FIAT/regge.py | 4 +- FIAT/restricted.py | 17 +-- FIAT/serendipity.py | 2 +- FIAT/tensor_product.py | 4 +- test/unit/test_fiat.py | 28 +++- test/unit/test_hct.py | 4 +- test/unit/test_macro.py | 21 ++- test/unit/test_stokes_complex.py | 144 ++++++++++++++++++ 55 files changed, 1090 insertions(+), 398 deletions(-) create mode 100644 FIAT/alfeld_sorokina.py create mode 100644 FIAT/arnold_qin.py create mode 100644 FIAT/bernardi_raugel.py create mode 100644 FIAT/christiansen_hu.py create mode 100644 FIAT/guzman_neilan.py create mode 100644 test/unit/test_stokes_complex.py diff --git a/FIAT/P0.py b/FIAT/P0.py index 1a88ee2d5..01d9a6472 100644 --- a/FIAT/P0.py +++ b/FIAT/P0.py @@ -40,7 +40,7 @@ def __init__(self, ref_el): entity_ids[dim][entity] = [entity] if dim == sd else [] entity_permutations[dim][entity] = perms - super(P0Dual, self).__init__(nodes, ref_el, entity_ids, entity_permutations) + super().__init__(nodes, ref_el, entity_ids, entity_permutations) class P0(finite_element.CiarletElement): @@ -49,4 +49,4 @@ def __init__(self, ref_el): dual = P0Dual(ref_el) degree = 0 formdegree = ref_el.get_spatial_dimension() # n-form - super(P0, self).__init__(poly_set, dual, degree, formdegree) + super().__init__(poly_set, dual, degree, formdegree) diff --git a/FIAT/Sminus.py b/FIAT/Sminus.py index 6d9d4eaa3..506cd23fd 100644 --- a/FIAT/Sminus.py +++ b/FIAT/Sminus.py @@ -497,7 +497,7 @@ def __init__(self, ref_el, degree): self.basis = {(0, 0): Array(Sminus_list)} else: self.basis = {(0, 0, 0): Array(Sminus_list)} - super(TrimmedSerendipityEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") + super().__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") class TrimmedSerendipityFace(TrimmedSerendipity): @@ -526,4 +526,4 @@ def __init__(self, ref_el, degree): Sminus_list = EL + FL Sminus_list = [[-a[1], a[0]] for a in Sminus_list] self.basis = {(0, 0): Array(Sminus_list)} - super(TrimmedSerendipityFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") diff --git a/FIAT/SminusCurl.py b/FIAT/SminusCurl.py index ce42a194f..3aaab57bf 100644 --- a/FIAT/SminusCurl.py +++ b/FIAT/SminusCurl.py @@ -88,11 +88,11 @@ def __init__(self, ref_el, degree, mapping): entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) - super(TrimmedSerendipity, self).__init__(ref_el=ref_el, - dual=None, - order=degree, - formdegree=formdegree, - mapping=mapping) + super().__init__(ref_el=ref_el, + dual=None, + order=degree, + formdegree=formdegree, + mapping=mapping) topology = ref_el.get_topology() unflattening_map = compute_unflattening_map(topology) @@ -231,7 +231,7 @@ def __init__(self, ref_el, degree): Sminus_list = EL + FL self.basis = {(0, 0): Array(Sminus_list)} - super(TrimmedSerendipityCurl, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") def e_lambda_1_3d(deg, dx, dy, dz, x_mid, y_mid, z_mid): diff --git a/FIAT/SminusDiv.py b/FIAT/SminusDiv.py index 70527a09b..3be143fa2 100644 --- a/FIAT/SminusDiv.py +++ b/FIAT/SminusDiv.py @@ -78,11 +78,11 @@ def __init__(self, ref_el, degree, mapping): entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) - super(TrimmedSerendipity, self).__init__(ref_el=ref_el, - dual=None, - order=degree, - formdegree=formdegree, - mapping=mapping) + super().__init__(ref_el=ref_el, + dual=None, + order=degree, + formdegree=formdegree, + mapping=mapping) topology = ref_el.get_topology() unflattening_map = compute_unflattening_map(topology) @@ -198,7 +198,7 @@ def __init__(self, ref_el, degree): IL = () Sminus_list = FL + IL self.basis = {(0, 0, 0): Array(Sminus_list)} - super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") else: # Put all 2 dimensional stuff here. @@ -213,7 +213,7 @@ def __init__(self, ref_el, degree): Sminus_list = EL + FL Sminus_list = [[-a[1], a[0]] for a in Sminus_list] self.basis = {(0, 0): Array(Sminus_list)} - super(TrimmedSerendipityDiv, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") def f_lambda_2_3d(degree, dx, dy, dz, x_mid, y_mid, z_mid): diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 965bf63aa..65afac645 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -7,9 +7,14 @@ # Import finite element classes from FIAT.finite_element import FiniteElement, CiarletElement # noqa: F401 from FIAT.argyris import Argyris +from FIAT.bernardi_raugel import BernardiRaugel from FIAT.bernstein import Bernstein from FIAT.bell import Bell from FIAT.hct import HsiehCloughTocher +from FIAT.alfeld_sorokina import AlfeldSorokina +from FIAT.arnold_qin import ArnoldQin +from FIAT.guzman_neilan import GuzmanNeilan +from FIAT.christiansen_hu import ChristiansenHu from FIAT.johnson_mercier import JohnsonMercier from FIAT.brezzi_douglas_marini import BrezziDouglasMarini from FIAT.Sminus import TrimmedSerendipityEdge # noqa: F401 @@ -31,8 +36,7 @@ from FIAT.morley import Morley from FIAT.nedelec import Nedelec from FIAT.nedelec_second_kind import NedelecSecondKind -from FIAT.powell_sabin import QuadraticPowellSabin6 -from FIAT.powell_sabin import QuadraticPowellSabin12 +from FIAT.powell_sabin import QuadraticPowellSabin6, QuadraticPowellSabin12 from FIAT.hierarchical import Legendre, IntegratedLegendre from FIAT.P0 import P0 from FIAT.raviart_thomas import RaviartThomas @@ -65,6 +69,7 @@ # List of supported elements and mapping to element classes supported_elements = {"Argyris": Argyris, "Bell": Bell, + "Bernardi-Raugel": BernardiRaugel, "Bernstein": Bernstein, "Brezzi-Douglas-Marini": BrezziDouglasMarini, "Brezzi-Douglas-Fortin-Marini": BrezziDouglasFortinMarini, @@ -79,7 +84,11 @@ "Discontinuous Taylor": DiscontinuousTaylor, "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas, "Hermite": CubicHermite, - "HsiehCloughTocher": HsiehCloughTocher, + "Hsieh-Clough-Tocher": HsiehCloughTocher, + "Alfeld-Sorokina": AlfeldSorokina, + "Arnold-Qin": ArnoldQin, + "Christiansen-Hu": ChristiansenHu, + "Guzman-Neilan": GuzmanNeilan, "Johnson-Mercier": JohnsonMercier, "Lagrange": Lagrange, "Kong-Mulder-Veldhuizen": KongMulderVeldhuizen, diff --git a/FIAT/alfeld_sorokina.py b/FIAT/alfeld_sorokina.py new file mode 100644 index 000000000..13ef25bce --- /dev/null +++ b/FIAT/alfeld_sorokina.py @@ -0,0 +1,106 @@ +# Copyright (C) 2024 Pablo D. Brubeck +# +# This file is part of FIAT (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# +# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024 + +from FIAT import finite_element, dual_set, polynomial_set +from FIAT.functional import ComponentPointEvaluation, PointDivergence, IntegralMoment +from FIAT.quadrature_schemes import create_quadrature +from FIAT.quadrature import FacetQuadratureRule +from FIAT.macro import CkPolynomialSet, AlfeldSplit + +import numpy + + +def AlfeldSorokinaSpace(ref_el, degree): + """Return a vector-valued C^0 PolynomialSet on an Alfeld split whose + divergence is also C^0. This works on any simplex and for all + polynomial degrees.""" + ref_complex = AlfeldSplit(ref_el) + sd = ref_complex.get_spatial_dimension() + C0 = CkPolynomialSet(ref_complex, degree, order=0, shape=(sd,), variant="bubble") + expansion_set = C0.get_expansion_set() + num_members = C0.get_num_members() + coeffs = C0.get_coeffs() + + facet_el = ref_complex.construct_subelement(sd-1) + phi = polynomial_set.ONPolynomialSet(facet_el, 0 if sd == 1 else degree-1) + Q = create_quadrature(facet_el, 2 * phi.degree) + qpts, qwts = Q.get_points(), Q.get_weights() + phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)] + weights = numpy.multiply(phi_at_qpts, qwts) + + rows = [] + for facet in ref_complex.get_interior_facets(sd-1): + n = ref_complex.compute_normal(facet) + jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=1) + div_jump = n[:, None, None] * jumps[1][None, ...] + r = numpy.tensordot(div_jump, weights, axes=(-1, -1)) + rows.append(r.reshape(num_members, -1).T) + + if len(rows) > 0: + dual_mat = numpy.vstack(rows) + _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) + tol = sig[0] * 1E-10 + num_sv = len([s for s in sig if abs(s) > tol]) + coeffs = numpy.tensordot(vt[num_sv:], coeffs, axes=(-1, 0)) + return polynomial_set.PolynomialSet(ref_complex, degree, degree, expansion_set, coeffs) + + +class AlfeldSorokinaDualSet(dual_set.DualSet): + def __init__(self, ref_el, degree, variant=None): + if degree != 2: + raise NotImplementedError("Alfeld-Sorokina only defined for degree = 2") + if variant is None: + variant = "integral" + if variant not in {"integral", "point"}: + raise ValueError(f"Invalid variant {variant}") + + top = ref_el.get_topology() + sd = ref_el.get_spatial_dimension() + entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} + + nodes = [] + dims = (0, 1) if variant == "point" else (0,) + for dim in dims: + for entity in sorted(top[dim]): + cur = len(nodes) + pts = ref_el.make_points(dim, entity, degree) + if dim == 0: + pt, = pts + nodes.append(PointDivergence(ref_el, pt)) + nodes.extend(ComponentPointEvaluation(ref_el, k, (sd,), pt) + for pt in pts for k in range(sd)) + entity_ids[dim][entity].extend(range(cur, len(nodes))) + + if variant == "integral": + # Edge moments against quadratic Legendre (mean-free bubble) + dim = 1 + facet = ref_el.construct_subelement(dim) + Q = create_quadrature(facet, degree+dim+1) + f_at_qpts = facet.compute_bubble(Q.get_points()) + f_at_qpts -= numpy.dot(f_at_qpts, Q.get_weights()) / facet.volume() + for entity in sorted(top[dim]): + cur = len(nodes) + Q_mapped = FacetQuadratureRule(ref_el, dim, entity, Q) + detJ = Q_mapped.jacobian_determinant() + phi = f_at_qpts / detJ + nodes.extend(IntegralMoment(ref_el, Q_mapped, phi, comp=k, shp=(sd,)) + for k in range(sd)) + entity_ids[dim][entity].extend(range(cur, len(nodes))) + + super().__init__(nodes, ref_el, entity_ids) + + +class AlfeldSorokina(finite_element.CiarletElement): + """The Alfeld-Sorokina C^0 quadratic macroelement with C^0 divergence. This element + belongs to a Stokes complex, and is paired with Lagrange(ref_el, 1, variant="alfeld").""" + def __init__(self, ref_el, degree=2): + dual = AlfeldSorokinaDualSet(ref_el, degree) + poly_set = AlfeldSorokinaSpace(ref_el, degree) + formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form + super().__init__(poly_set, dual, degree, formdegree, + mapping="contravariant piola") diff --git a/FIAT/argyris.py b/FIAT/argyris.py index ba35b61d7..d98dbbe3b 100644 --- a/FIAT/argyris.py +++ b/FIAT/argyris.py @@ -82,7 +82,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): entity_ids[2][0] = list(range(cur, len(nodes))) else: raise ValueError("Invalid variant for Argyris") - super(ArgyrisDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class Argyris(finite_element.CiarletElement): @@ -107,4 +107,4 @@ def __init__(self, ref_el, degree=5, variant=None): poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant="bubble") dual = ArgyrisDualSet(ref_el, degree, variant, interpolant_deg) - super(Argyris, self).__init__(poly_set, dual, degree) + super().__init__(poly_set, dual, degree) diff --git a/FIAT/arnold_qin.py b/FIAT/arnold_qin.py new file mode 100644 index 000000000..c2a421dd5 --- /dev/null +++ b/FIAT/arnold_qin.py @@ -0,0 +1,71 @@ +# This file is part of FIAT (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# +# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024 + +from FIAT import finite_element, polynomial_set +from FIAT.bernardi_raugel import BernardiRaugelDualSet +from FIAT.quadrature_schemes import create_quadrature +from FIAT.reference_element import TRIANGLE +from FIAT.macro import CkPolynomialSet +from FIAT.hct import HsiehCloughTocher + +import numpy + + +def ArnoldQinSpace(ref_el, degree, reduced=False): + """Return a basis for the Arnold-Qin space + curl(HCT-red) + P_0 x if reduced = True, and + curl(HCT) + P_0 x if reduced = False.""" + if ref_el.get_shape() != TRIANGLE: + raise ValueError("Arnold-Qin only defined on triangles") + if degree != 2: + raise ValueError("Arnold-Qin only defined for degree = 2") + sd = ref_el.get_spatial_dimension() + HCT = HsiehCloughTocher(ref_el, degree+1, reduced=True) + ref_complex = HCT.get_reference_complex() + Q = create_quadrature(ref_complex, 2 * degree) + qpts, qwts = Q.get_points(), Q.get_weights() + + x = qpts.T + bary = numpy.asarray(ref_el.make_points(sd, 0, sd+1)) + P0x_at_qpts = x[None, :, :] - bary[:, :, None] + + tab = HCT.tabulate(1, qpts) + curl_at_qpts = numpy.stack([tab[(0, 1)], -tab[(1, 0)]], axis=1) + if reduced: + curl_at_qpts = curl_at_qpts[:9] + + C0 = CkPolynomialSet(ref_complex, degree, order=0, scale=1, variant="bubble") + C0_at_qpts = C0.tabulate(qpts)[(0,) * sd] + duals = numpy.multiply(C0_at_qpts, qwts) + M = numpy.dot(duals, C0_at_qpts.T) + duals = numpy.linalg.solve(M, duals) + + # Remove the constant nullspace + ids = [0, 3, 6] + A = numpy.asarray([[1, 1, 1], [1, -1, 0], [0, -1, 1]]) + phis = curl_at_qpts + phis[ids] = numpy.tensordot(A, phis[ids], axes=(-1, 0)) + # Replace the constant nullspace with P_0 x + phis[0] = P0x_at_qpts + coeffs = numpy.tensordot(phis, duals, axes=(-1, -1)) + return polynomial_set.PolynomialSet(ref_complex, degree, degree, + C0.get_expansion_set(), coeffs) + + +class ArnoldQin(finite_element.CiarletElement): + """The Arnold-Qin C^0(Alfeld) quadratic macroelement with divergence in P0. + This element belongs to a Stokes complex, and is paired with unsplit DG0.""" + def __init__(self, ref_el, degree=2, reduced=False): + poly_set = ArnoldQinSpace(ref_el, degree) + if reduced: + subdegree = 1 + mapping = "contravariant piola" + else: + subdegree = degree + mapping = "affine" + dual = BernardiRaugelDualSet(ref_el, degree, subdegree) + formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form + super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index 5e8bdbbd1..5c12fcc76 100644 --- a/FIAT/arnold_winther.py +++ b/FIAT/arnold_winther.py @@ -25,7 +25,7 @@ class ArnoldWintherNCDual(DualSet): - def __init__(self, cell, degree): + def __init__(self, cell, degree=2): if not degree == 2: raise ValueError("Nonconforming Arnold-Winther elements are" "only defined for degree 2.") @@ -70,24 +70,23 @@ def __init__(self, cell, degree): dof_ids[1][entity_id].append(dof_cur) dof_cur += 1 - super(ArnoldWintherNCDual, self).__init__(dofs, cell, dof_ids) + super().__init__(dofs, cell, dof_ids) class ArnoldWintherNC(CiarletElement): """The definition of the nonconforming Arnold-Winther element. """ - def __init__(self, cell, degree): + def __init__(self, cell, degree=2): assert degree == 2, "Only defined for degree 2" Ps = ONSymTensorPolynomialSet(cell, degree) Ls = ArnoldWintherNCDual(cell, degree) mapping = "double contravariant piola" - super(ArnoldWintherNC, self).__init__(Ps, Ls, degree, - mapping=mapping) + super().__init__(Ps, Ls, degree, mapping=mapping) class ArnoldWintherDual(DualSet): - def __init__(self, cell, degree): + def __init__(self, cell, degree=3): if not degree == 3: raise ValueError("Arnold-Winther elements are" "only defined for degree 3.") @@ -152,15 +151,15 @@ def __init__(self, cell, degree): dof_ids[2][0] += list(range(dof_cur, dof_cur+6)) - super(ArnoldWintherDual, self).__init__(dofs, cell, dof_ids) + super().__init__(dofs, cell, dof_ids) class ArnoldWinther(CiarletElement): """The definition of the conforming Arnold-Winther element. """ - def __init__(self, cell, degree): + def __init__(self, cell, degree=3): assert degree == 3, "Only defined for degree 3" Ps = ONSymTensorPolynomialSet(cell, degree) Ls = ArnoldWintherDual(cell, degree) mapping = "double contravariant piola" - super(ArnoldWinther, self).__init__(Ps, Ls, degree, mapping=mapping) + super().__init__(Ps, Ls, degree, mapping=mapping) diff --git a/FIAT/barycentric_interpolation.py b/FIAT/barycentric_interpolation.py index 7e2f4111d..25862bf2f 100644 --- a/FIAT/barycentric_interpolation.py +++ b/FIAT/barycentric_interpolation.py @@ -73,7 +73,7 @@ def __init__(self, ref_el, pts): self.degree = max(len(wts) for wts in self.weights)-1 self.recurrence_order = self.degree + 1 - super(LagrangeLineExpansionSet, self).__init__(ref_el) + super().__init__(ref_el) self.continuity = None if len(self.x) == sum(len(xk) for xk in self.nodes) else "C0" def get_num_members(self, n): @@ -123,5 +123,4 @@ def __init__(self, ref_el, pts, shape=tuple()): coeffs[cur_idx] = 1.0 cur_bf += 1 - super(LagrangePolynomialSet, self).__init__(ref_el, degree, embedded_degree, - expansion_set, coeffs) + super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs) diff --git a/FIAT/bell.py b/FIAT/bell.py index 2603d3a77..1d3fe11b3 100644 --- a/FIAT/bell.py +++ b/FIAT/bell.py @@ -53,7 +53,7 @@ def __init__(self, ref_el): from FIAT.jacobi import eval_jacobi rline = ufc_simplex(1) q1d = create_quadrature(rline, 8) - q1dpts = q1d.get_points() + q1dpts = q1d.get_points()[:, 0] leg4_at_qpts = eval_jacobi(0, 0, 4, 2.0*q1dpts - 1) imond = functional.IntegralMomentOfNormalDerivative @@ -64,7 +64,7 @@ def __init__(self, ref_el): entity_ids[2] = {0: []} - super(BellDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class Bell(finite_element.CiarletElement): @@ -73,4 +73,4 @@ class Bell(finite_element.CiarletElement): def __init__(self, ref_el): poly_set = polynomial_set.ONPolynomialSet(ref_el, 5) dual = BellDualSet(ref_el) - super(Bell, self).__init__(poly_set, dual, 5) + super().__init__(poly_set, dual, 5) diff --git a/FIAT/bernardi_raugel.py b/FIAT/bernardi_raugel.py new file mode 100644 index 000000000..68e8720f2 --- /dev/null +++ b/FIAT/bernardi_raugel.py @@ -0,0 +1,102 @@ +# This file is part of FIAT (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# +# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024 + +# This is not quite Bernardi-Raugel, but it has 2*dim*(dim+1) dofs and includes +# dim**2-1 extra constraint functionals. The first (dim+1)**2 basis functions +# are the reference element bfs, but the extra dim**2-1 are used in the +# transformation theory. + +from FIAT import finite_element, dual_set, polynomial_set, expansions +from FIAT.functional import ComponentPointEvaluation, FrobeniusIntegralMoment +from FIAT.quadrature_schemes import create_quadrature +from FIAT.quadrature import FacetQuadratureRule +import numpy + + +def ExtendedBernardiRaugelSpace(ref_el, subdegree): + r"""Return a basis for the extended Bernardi-Raugel space. + P_k^d + (P_{dim} \ P_{dim-1})^d""" + sd = ref_el.get_spatial_dimension() + if subdegree >= sd: + raise ValueError("The Bernardi-Raugel space is only defined for subdegree < dim") + Pd = polynomial_set.ONPolynomialSet(ref_el, sd, shape=(sd,), scale=1, variant="bubble") + dimPd = expansions.polynomial_dimension(ref_el, sd, continuity="C0") + entity_ids = expansions.polynomial_entity_ids(ref_el, sd, continuity="C0") + ids = [i + j * dimPd + for dim in (*tuple(range(subdegree)), sd-1) + for f in sorted(entity_ids[dim]) + for i in entity_ids[dim][f] + for j in range(sd)] + return Pd.take(ids) + + +class BernardiRaugelDualSet(dual_set.DualSet): + """The Bernardi-Raugel dual set.""" + def __init__(self, ref_complex, degree, subdegree=1, reduced=False): + ref_el = ref_complex.get_parent() or ref_complex + sd = ref_el.get_spatial_dimension() + if subdegree > sd: + raise ValueError("The Bernardi-Raugel dual is only defined for subdegree <= dim") + top = ref_el.get_topology() + entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} + + # Point evaluation at lattice points + nodes = [] + for dim in range(subdegree): + for entity in sorted(top[dim]): + cur = len(nodes) + pts = ref_el.make_points(dim, entity, subdegree) + nodes.extend(ComponentPointEvaluation(ref_el, comp, (sd,), pt) + for pt in pts for comp in range(sd)) + entity_ids[dim][entity].extend(range(cur, len(nodes))) + + if subdegree < sd: + # Face moments of normal/tangential components against mean-free bubbles + facet = ref_complex.construct_subcomplex(sd-1) + Q = create_quadrature(facet, 2*degree) + if degree == 1 and facet.is_macrocell(): + P = polynomial_set.ONPolynomialSet(facet, degree, scale=1, variant="bubble") + f_at_qpts = P.tabulate(Q.get_points())[(0,)*(sd-1)][-1] + else: + ref_facet = facet.get_parent() or facet + f_at_qpts = ref_facet.compute_bubble(Q.get_points()) + f_at_qpts -= numpy.dot(f_at_qpts, Q.get_weights()) / facet.volume() + + Qs = {f: FacetQuadratureRule(ref_el, sd-1, f, Q) + for f in sorted(top[sd-1])} + + thats = {f: ref_el.compute_tangents(sd-1, f) + for f in sorted(top[sd-1])} + + R = numpy.array([[0, 1], [-1, 0]]) + ndir = 1 if reduced else sd + for i in range(ndir): + for f, Q_mapped in Qs.items(): + cur = len(nodes) + if i == 0: + udir = numpy.dot(R, *thats[f]) if sd == 2 else numpy.cross(*thats[f]) + else: + udir = thats[f][i-1] + detJ = Q_mapped.jacobian_determinant() + phi_at_qpts = udir[:, None] * f_at_qpts[None, :] / detJ + nodes.append(FrobeniusIntegralMoment(ref_el, Q_mapped, phi_at_qpts)) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) + + super().__init__(nodes, ref_el, entity_ids) + + +class BernardiRaugel(finite_element.CiarletElement): + """The Bernardi-Raugel extended element.""" + def __init__(self, ref_el, degree=None, subdegree=1): + sd = ref_el.get_spatial_dimension() + if degree is None: + degree = sd + if degree != sd: + raise ValueError("Bernardi-Raugel only defined for degree = dim") + poly_set = ExtendedBernardiRaugelSpace(ref_el, subdegree) + dual = BernardiRaugelDualSet(ref_el, degree, subdegree=subdegree) + formdegree = sd - 1 # (n-1)-form + super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index 81768d12d..096e10f16 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -43,7 +43,7 @@ def __init__(self, ref_el, degree): # Leave nodes unimplemented for now nodes.append(None) - super(BernsteinDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class Bernstein(FiniteElement): @@ -52,7 +52,7 @@ class Bernstein(FiniteElement): def __init__(self, ref_el, degree): dual = BernsteinDualSet(ref_el, degree) k = 0 # 0-form - super(Bernstein, self).__init__(ref_el, dual, degree, k) + super().__init__(ref_el, dual, degree, k) def degree(self): """The degree of the polynomial space.""" diff --git a/FIAT/brezzi_douglas_fortin_marini.py b/FIAT/brezzi_douglas_fortin_marini.py index fb8f81bd8..4361e70f6 100644 --- a/FIAT/brezzi_douglas_fortin_marini.py +++ b/FIAT/brezzi_douglas_fortin_marini.py @@ -58,7 +58,7 @@ def __init__(self, ref_el, degree): entity_ids[sd] = {0: list(range(cur, cur + tangent_count))} cur += tangent_count - super(BDFMDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) def BDFMSpace(ref_el, order): @@ -112,5 +112,5 @@ def __init__(self, ref_el, degree): poly_set = BDFMSpace(ref_el, degree) dual = BDFMDualSet(ref_el, degree - 1) formdegree = ref_el.get_spatial_dimension() - 1 - super(BrezziDouglasFortinMarini, self).__init__(poly_set, dual, degree, formdegree, - mapping="contravariant piola") + super().__init__(poly_set, dual, degree, formdegree, + mapping="contravariant piola") diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 2c9aeb892..0014b5277 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -65,7 +65,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): for phi in Ned_at_qpts) entity_ids[sd][0] = list(range(cur, len(nodes))) - super(BDMDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class BrezziDouglasMarini(finite_element.CiarletElement): @@ -99,5 +99,4 @@ def __init__(self, ref_el, degree, variant=None): poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd, )) dual = BDMDualSet(ref_el, degree, variant, interpolant_deg) formdegree = sd - 1 # (n-1)-form - super(BrezziDouglasMarini, self).__init__(poly_set, dual, degree, formdegree, - mapping="contravariant piola") + super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 7a5ab500c..786a5d309 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -76,9 +76,9 @@ def __init__(self, ref_el, degree, mapping): entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) # Set up FiniteElement - super(BrezziDouglasMariniCube, self).__init__(ref_el=ref_el, dual=None, - order=degree, formdegree=1, - mapping=mapping) + super().__init__(ref_el=ref_el, dual=None, + order=degree, formdegree=1, + mapping=mapping) # Store unflattened entity ID dictionaries topology = ref_el.get_topology() @@ -272,8 +272,7 @@ def __init__(self, ref_el, degree): bdmce_list = construct_bdmce_basis(ref_el, degree) self.basis = {(0, 0): Array(bdmce_list)} - super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, degree=degree, - mapping="covariant piola") + super().__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") class BrezziDouglasMariniCubeFace(BrezziDouglasMariniCube): @@ -292,8 +291,7 @@ def __init__(self, ref_el, degree): bdmcf_list = [[-a[1], a[0]] for a in bdmce_list] self.basis = {(0, 0): Array(bdmcf_list)} - super(BrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree, - mapping="contravariant piola") + super().__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") def construct_bdmce_basis(ref_el, degree): diff --git a/FIAT/bubble.py b/FIAT/bubble.py index 1747d9e5c..0de7b2566 100644 --- a/FIAT/bubble.py +++ b/FIAT/bubble.py @@ -23,18 +23,18 @@ def __init__(self, ref_el, degree, codim): if len(dofs) == 0: raise RuntimeError('Bubble element of degree %d and codimension %d has no dofs' % (degree, codim)) - super(CodimBubble, self).__init__(element, indices=dofs) + super().__init__(element, indices=dofs) 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(Bubble, self).__init__(ref_el, degree, codim=0) + super().__init__(ref_el, degree, codim=0) 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(FacetBubble, self).__init__(ref_el, degree, codim=1) + super().__init__(ref_el, degree, codim=1) diff --git a/FIAT/check_format_variant.py b/FIAT/check_format_variant.py index b42c01ac6..1431876de 100644 --- a/FIAT/check_format_variant.py +++ b/FIAT/check_format_variant.py @@ -1,6 +1,6 @@ import re -from FIAT.macro import AlfeldSplit, IsoSplit +from FIAT.macro import IsoSplit, AlfeldSplit, WorseyFarinSplit, PowellSabinSplit, PowellSabin12Split # dicts mapping Lagrange variant names to recursivenodes family names supported_cg_variants = { @@ -17,6 +17,14 @@ "gll": "gll", "gl": "gl"} +supported_splits = { + "iso": IsoSplit, + "alfeld": AlfeldSplit, + "worsey-farin": WorseyFarinSplit, + "powell-sabin": PowellSabinSplit, + "powell-sabin(12)": PowellSabin12Split, +} + def check_format_variant(variant, degree): if variant is None: @@ -44,7 +52,7 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False): variant may be a single option or comma-separated pair indicating the dof type (integral, equispaced, spectral, etc) - and the type of splitting to give a macro-element (Alfeld, iso) + and the type of splitting to give a macro-element (Alfeld, Powell-Sabin, iso) """ if variant is None: variant = "integral" if integral else "equispaced" @@ -66,10 +74,8 @@ def parse_lagrange_variant(variant, discontinuous=False, integral=False): for pre_opt in options: opt = pre_opt.lower() - if opt == "alfeld": - splitting = AlfeldSplit - elif opt == "iso": - splitting = IsoSplit + if opt in supported_splits: + splitting = supported_splits[opt] elif opt.startswith("iso"): match = re.match(r"^iso(?:\((\d+)\))?$", opt) k, = match.groups() diff --git a/FIAT/christiansen_hu.py b/FIAT/christiansen_hu.py new file mode 100644 index 000000000..d9d5b33e0 --- /dev/null +++ b/FIAT/christiansen_hu.py @@ -0,0 +1,78 @@ +# This file is part of FIAT (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# +# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024 + +# This is not quite Christiansen-Hu, but it has 2*dim*(dim+1) dofs and includes +# dim**2-1 extra constraint functionals. The first (dim+1)**2 basis functions +# are the reference element bfs, but the extra dim**2-1 are used in the +# transformation theory. + +from FIAT import finite_element, polynomial_set +from FIAT.bernardi_raugel import BernardiRaugelDualSet +from FIAT.quadrature_schemes import create_quadrature +from FIAT.macro import CkPolynomialSet, WorseyFarinSplit + +import numpy + + +def ChristiansenHuSpace(ref_el, degree, reduced=False): + """Return a basis for the Christianse-Hu space + set(v in C0 P1(WF)^d : div(v) = 0) + P_0 x if reduced = True, and + this space is agumented with rotated facet bubbles if reduced = False.""" + sd = ref_el.get_spatial_dimension() + ref_complex = WorseyFarinSplit(ref_el) + C0 = CkPolynomialSet(ref_complex, degree, order=0, shape=(sd,), scale=1, variant="bubble") + Q = create_quadrature(ref_complex, degree-1) + tab = C0.tabulate(Q.get_points(), 1) + divC0 = sum(tab[alpha][:, alpha.index(1), :] for alpha in tab if sum(alpha) == 1) + + _, sig, vt = numpy.linalg.svd(divC0.T, full_matrices=True) + tol = sig[0] * 1E-10 + num_sv = len([s for s in sig if abs(s) > tol]) + coeffs = numpy.tensordot(vt[num_sv:], C0.get_coeffs(), axes=(-1, 0)) + + verts = numpy.array(ref_complex.get_vertices()) + WT = verts[-1] + P0x_coeffs = numpy.transpose(verts - WT[None, :]) + coeffs = numpy.concatenate((coeffs, P0x_coeffs[None, ...]), axis=0) + + if not reduced: + # Compute the primal basis via Vandermonde and extract the facet bubbles + dual = BernardiRaugelDualSet(ref_complex, degree, reduced=True) + dualmat = dual.to_riesz(C0) + V = numpy.tensordot(dualmat, coeffs, axes=((1, 2), (1, 2))) + coeffs = numpy.tensordot(numpy.linalg.inv(V.T), coeffs, axes=(-1, 0)) + facet_bubbles = coeffs[-(sd+1):] + + # Rotate the facet bubbles onto the tangent space of the facet + # NOTE they are not aligned with the normal, but they point in the direction + # that connects the split point on the facet with the split point of the cell + WF = verts[sd+1:-1] + top = ref_el.get_topology() + ext = [] + for f in top[sd-1]: + ehat = WF[f] - WT + FB = numpy.dot(ehat, facet_bubbles[f]) + thats = ref_el.compute_tangents(sd-1, f) + for that in thats: + ext.append(that[:, None] * FB[None, :]) + ext_coeffs = numpy.array(ext) + coeffs = numpy.concatenate((coeffs, ext_coeffs), axis=0) + + return polynomial_set.PolynomialSet(ref_complex, degree, degree, + C0.get_expansion_set(), coeffs) + + +class ChristiansenHu(finite_element.CiarletElement): + """The Christiansen-Hu C^0(Worsey-Farin) linear macroelement with divergence in P0. + This element belongs to a Stokes complex, and is paired with unsplit DG0.""" + def __init__(self, ref_el, degree=1): + if degree != 1: + raise ValueError("Christiansen-Hu only defined for degree = 1") + poly_set = ChristiansenHuSpace(ref_el, degree) + ref_complex = poly_set.get_reference_element() + dual = BernardiRaugelDualSet(ref_complex, degree) + formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form + super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/crouzeix_raviart.py b/FIAT/crouzeix_raviart.py index b8fa3f80e..b809f7a87 100644 --- a/FIAT/crouzeix_raviart.py +++ b/FIAT/crouzeix_raviart.py @@ -46,7 +46,7 @@ def __init__(self, cell, degree): entity_ids[d - 1][i] += [i] # Initialize super-class - super(CrouzeixRaviartDualSet, self).__init__(nodes, cell, entity_ids) + super().__init__(nodes, cell, entity_ids) class CrouzeixRaviart(finite_element.CiarletElement): @@ -67,4 +67,4 @@ def __init__(self, cell, degree): # FiniteElement space = polynomial_set.ONPolynomialSet(cell, 1) dual = CrouzeixRaviartDualSet(cell, 1) - super(CrouzeixRaviart, self).__init__(space, dual, 1) + super().__init__(space, dual, 1) diff --git a/FIAT/discontinuous_lagrange.py b/FIAT/discontinuous_lagrange.py index e086421e5..fd08f8fa3 100644 --- a/FIAT/discontinuous_lagrange.py +++ b/FIAT/discontinuous_lagrange.py @@ -169,7 +169,7 @@ def __init__(self, ref_el, degree, point_variant="equispaced"): entity_permutations[dim][entity] = perms entity_ids[dim][0] = list(range(len(nodes))) - super(BrokenLagrangeDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) + super().__init__(nodes, ref_el, entity_ids, entity_permutations) class DiscontinuousLagrangeDualSet(dual_set.DualSet): @@ -196,7 +196,7 @@ def __init__(self, ref_el, degree, point_variant="equispaced"): degree, variant=point_variant) nodes.extend(functional.PointEvaluation(ref_el, x) for x in pts) entity_ids[dim][entity] = list(range(cur, len(nodes))) - super(DiscontinuousLagrangeDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) + super().__init__(nodes, ref_el, entity_ids, entity_permutations) class DiscontinuousLagrange(finite_element.CiarletElement): @@ -220,7 +220,7 @@ def __new__(cls, ref_el, degree, variant="equispaced"): 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(DiscontinuousLagrange, cls).__new__(cls) + return super().__new__(cls) def __init__(self, ref_el, degree, variant="equispaced"): splitting, point_variant = parse_lagrange_variant(variant, discontinuous=True) @@ -238,4 +238,4 @@ def __init__(self, ref_el, degree, variant="equispaced"): else: poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) formdegree = ref_el.get_spatial_dimension() # n-form - super(DiscontinuousLagrange, self).__init__(poly_set, dual, degree, formdegree) + super().__init__(poly_set, dual, degree, formdegree) diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index af02f8861..dbb566906 100644 --- a/FIAT/discontinuous_pc.py +++ b/FIAT/discontinuous_pc.py @@ -39,11 +39,11 @@ def __init__(self, ref_el): dual.entity_permutations = None degree = 0 formdegree = ref_el.get_spatial_dimension() # n-form - super(DPC0, self).__init__(poly_set=poly_set, - dual=dual, - order=degree, - ref_complex=ref_el, - formdegree=formdegree) + super().__init__(poly_set=poly_set, + dual=dual, + order=degree, + ref_complex=ref_el, + formdegree=formdegree) class DPCDualSet(dual_set.DualSet): @@ -94,7 +94,7 @@ def __init__(self, ref_el, flat_el, degree): entity_ids[dim][entity] = [] entity_ids[dim][0] = list(range(len(nodes))) - super(DPCDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class HigherOrderDPC(finite_element.CiarletElement): @@ -105,11 +105,11 @@ def __init__(self, ref_el, degree): poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[flat_el], degree) dual = DPCDualSet(ref_el, flat_el, degree) formdegree = flat_el.get_spatial_dimension() # n-form - super(HigherOrderDPC, self).__init__(poly_set=poly_set, - dual=dual, - order=degree, - ref_complex=ref_el, - formdegree=formdegree) + super().__init__(poly_set=poly_set, + dual=dual, + order=degree, + ref_complex=ref_el, + formdegree=formdegree) def DPC(ref_el, degree): diff --git a/FIAT/discontinuous_raviart_thomas.py b/FIAT/discontinuous_raviart_thomas.py index 12a73e650..8e3444e94 100644 --- a/FIAT/discontinuous_raviart_thomas.py +++ b/FIAT/discontinuous_raviart_thomas.py @@ -49,7 +49,7 @@ def __init__(self, ref_el, degree): # cell dofs entity_ids[sd] = {0: list(range(len(nodes)))} - super(DRTDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class DiscontinuousRaviartThomas(finite_element.CiarletElement): @@ -59,5 +59,4 @@ def __init__(self, ref_el, degree): poly_set = RTSpace(ref_el, degree) dual = DRTDualSet(ref_el, degree) - super(DiscontinuousRaviartThomas, self).__init__(poly_set, dual, degree, - mapping="contravariant piola") + super().__init__(poly_set, dual, degree, mapping="contravariant piola") diff --git a/FIAT/discontinuous_taylor.py b/FIAT/discontinuous_taylor.py index 7c27fd596..1438d1d3b 100644 --- a/FIAT/discontinuous_taylor.py +++ b/FIAT/discontinuous_taylor.py @@ -36,7 +36,7 @@ def __init__(self, ref_el, degree): for d in range(dim + 1)} entity_ids[dim][0] = list(range(len(nodes))) - super(DiscontinuousTaylorDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class HigherOrderDiscontinuousTaylor(finite_element.CiarletElement): @@ -46,7 +46,7 @@ def __init__(self, ref_el, degree): poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) dual = DiscontinuousTaylorDualSet(ref_el, degree) formdegree = ref_el.get_spatial_dimension() # n-form - super(HigherOrderDiscontinuousTaylor, self).__init__(poly_set, dual, degree, formdegree) + super().__init__(poly_set, dual, degree, formdegree) def DiscontinuousTaylor(ref_el, degree): diff --git a/FIAT/enriched.py b/FIAT/enriched.py index 97ac5c57c..b57da8401 100644 --- a/FIAT/enriched.py +++ b/FIAT/enriched.py @@ -59,7 +59,7 @@ def __init__(self, *elements): nodes = list(chain.from_iterable(e.dual_basis() for e in elements)) dual = DualSet(nodes, ref_el, entity_ids) - super(EnrichedElement, self).__init__(ref_el, dual, order, formdegree, mapping) + super().__init__(ref_el, dual, order, formdegree, mapping) # required degree (for quadrature) is definitely max self.polydegree = max(e.degree() for e in elements) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 8ff65f164..ac0eea173 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -247,7 +247,7 @@ def __new__(cls, *args, **kwargs): """Returns an ExpansionSet instance appropriate for the given reference element.""" if cls is not ExpansionSet: - return super(ExpansionSet, cls).__new__(cls) + return super().__new__(cls) try: ref_el = args[0] expansion_set = { @@ -518,7 +518,7 @@ class PointExpansionSet(ExpansionSet): def __init__(self, ref_el, **kwargs): if ref_el.get_spatial_dimension() != 0: raise ValueError("Must have a point") - super(PointExpansionSet, self).__init__(ref_el, **kwargs) + super().__init__(ref_el, **kwargs) def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): """Returns a dict of tabulations such that @@ -532,13 +532,13 @@ class LineExpansionSet(ExpansionSet): def __init__(self, ref_el, **kwargs): if ref_el.get_spatial_dimension() != 1: raise Exception("Must have a line") - super(LineExpansionSet, self).__init__(ref_el, **kwargs) + super().__init__(ref_el, **kwargs) def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): """Returns a dict of tabulations such that tabulations[alpha][i, j] = D^alpha phi_i(pts[j]).""" if self.variant is not None: - return super(LineExpansionSet, self)._tabulate_on_cell(n, pts, order=order, cell=cell, direction=direction) + return super()._tabulate_on_cell(n, pts, order=order, cell=cell, direction=direction) A, b = self.affine_mappings[cell] Jinv = A[0, 0] if direction is None else numpy.dot(A, direction) @@ -562,7 +562,7 @@ class TriangleExpansionSet(ExpansionSet): def __init__(self, ref_el, **kwargs): if ref_el.get_spatial_dimension() != 2: raise Exception("Must have a triangle") - super(TriangleExpansionSet, self).__init__(ref_el, **kwargs) + super().__init__(ref_el, **kwargs) class TetrahedronExpansionSet(ExpansionSet): @@ -570,7 +570,7 @@ class TetrahedronExpansionSet(ExpansionSet): def __init__(self, ref_el, **kwargs): if ref_el.get_spatial_dimension() != 3: raise Exception("Must be a tetrahedron") - super(TetrahedronExpansionSet, self).__init__(ref_el, **kwargs) + super().__init__(ref_el, **kwargs) def polynomial_dimension(ref_el, n, continuity=None): diff --git a/FIAT/fdm_element.py b/FIAT/fdm_element.py index 5f6408105..3c20c653b 100644 --- a/FIAT/fdm_element.py +++ b/FIAT/fdm_element.py @@ -140,7 +140,7 @@ def __init__(self, ref_el, degree, bc_order=1, formdegree=0, orthogonalize=False entity_permutations = {} entity_permutations[0] = {0: {0: []}, 1: {0: []}} entity_permutations[1] = {0: make_entity_permutations_simplex(1, degree + 1)} - super(FDMDual, self).__init__(nodes, ref_el, entity_ids, entity_permutations) + super().__init__(nodes, ref_el, entity_ids, entity_permutations) class FDMFiniteElement(finite_element.CiarletElement): @@ -171,7 +171,7 @@ def __init__(self, ref_el, degree): else: lr = quadrature.GaussLegendreQuadratureLineRule(ref_el, degree+1) poly_set = LagrangePolynomialSet(ref_el, lr.get_points()) - super(FDMFiniteElement, self).__init__(poly_set, dual, degree, self._formdegree) + super().__init__(poly_set, dual, degree, self._formdegree) class FDMLagrange(FDMFiniteElement): diff --git a/FIAT/functional.py b/FIAT/functional.py index a55fdb864..36989e607 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -12,13 +12,12 @@ # - a reference element domain # - type information -from collections import OrderedDict from itertools import chain import numpy import sympy -from FIAT import polynomial_set -from FIAT.quadrature import GaussLegendreQuadratureLineRule, QuadratureRule +from FIAT import polynomial_set, jacobi +from FIAT.quadrature import GaussLegendreQuadratureLineRule from FIAT.reference_element import UFCInterval as interval @@ -269,6 +268,16 @@ def __init__(self, ref_el, facet_no, pt): Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv") +class PointDivergence(Functional): + """Class representing point divergence of vector + functions at a particular point x.""" + + def __init__(self, ref_el, x): + dpt_dict = {x: [(1.0, alpha, (alpha.index(1),)) for alpha in polynomial_set.mis(len(x), 1)]} + + Functional.__init__(self, ref_el, (len(x),), {}, dpt_dict, "PointDiv") + + class IntegralMoment(Functional): """Functional representing integral of the input against some tabulated function f. @@ -313,17 +322,14 @@ def __init__(self, ref_el, facet_no, Q, f_at_qpts): sd = ref_el.get_spatial_dimension() # map points onto facet + transform = ref_el.get_entity_transform(sd-1, facet_no) + points = transform(Q.get_points()) + self.dpts = points + weights = numpy.multiply(f_at_qpts, Q.get_weights()) - fmap = ref_el.get_entity_transform(sd-1, facet_no) - qpts, qwts = Q.get_points(), Q.get_weights() - dpts = [fmap(pt) for pt in qpts] - self.dpts = dpts - - dpt_dict = OrderedDict() - - alphas = [tuple(1 if j == i else 0 for j in range(sd)) for i in range(sd)] - for j, pt in enumerate(dpts): - dpt_dict[tuple(pt)] = [(qwts[j]*n[i]*f_at_qpts[j], alphas[i], tuple()) for i in range(sd)] + alphas = tuple(map(tuple, numpy.eye(sd, dtype=int))) + dpt_dict = {tuple(pt): [(wt*n[i], alphas[i], tuple()) for i in range(sd)] + for pt, wt in zip(points, weights)} Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfNormalDerivative") @@ -337,20 +343,13 @@ def __init__(self, cell, s, entity, mom_deg, comp_deg, nm=""): shp = (sd,) quadpoints = comp_deg + 1 Q = GaussLegendreQuadratureLineRule(interval(), quadpoints) - legendre = numpy.polynomial.legendre.legval(2*Q.get_points()-1, [0]*mom_deg + [1]) - f_at_qpts = numpy.array([s*legendre[i] for i in range(quadpoints)]) - fmap = cell.get_entity_transform(sd-1, entity) - mappedqpts = [fmap(pt) for pt in Q.get_points()] - mappedQ = QuadratureRule(cell, mappedqpts, Q.get_weights()) - qwts = mappedQ.wts - qpts = mappedQ.pts - - pt_dict = OrderedDict() - - for k in range(len(qpts)): - pt_cur = tuple(qpts[k]) - pt_dict[pt_cur] = [(qwts[k] * f_at_qpts[k, i], (i,)) - for i in range(2)] + x = 2*Q.get_points()[:, 0]-1 + f_at_qpts = jacobi.eval_jacobi(0, 0, mom_deg, x) + transform = cell.get_entity_transform(sd-1, entity) + points = transform(Q.get_points()) + weights = numpy.multiply(f_at_qpts, Q.get_weights()) + pt_dict = {tuple(pt): [(wt*s[i], (i,)) for i in range(sd)] + for pt, wt in zip(points, weights)} super().__init__(cell, shp, pt_dict, {}, nm) @@ -377,32 +376,24 @@ def __init__(self, cell, s1, s2, entity, mom_deg, comp_deg, nm=""): # mom_deg is degree of moment, comp_deg is the total degree of # polynomial you might need to integrate (or something like that) sd = cell.get_spatial_dimension() - shp = (sd, sd) s1s2T = numpy.outer(s1, s2) + shp = s1s2T.shape quadpoints = comp_deg + 1 Q = GaussLegendreQuadratureLineRule(interval(), quadpoints) # The volume squared gets the Jacobian mapping from line interval # and the edge length into the functional. - legendre = numpy.polynomial.legendre.legval(2*Q.get_points()-1, [0]*mom_deg + [1]) * numpy.abs(cell.volume_of_subcomplex(1, entity))**2 - - f_at_qpts = numpy.array([s1s2T*legendre[i] for i in range(quadpoints)]) + x = 2*Q.get_points()[:, 0]-1 + f_at_qpts = jacobi.eval_jacobi(0, 0, mom_deg, x) * numpy.abs(cell.volume_of_subcomplex(1, entity))**2 # Map the quadrature points - fmap = cell.get_entity_transform(sd-1, entity) - mappedqpts = [fmap(pt) for pt in Q.get_points()] - mappedQ = QuadratureRule(cell, mappedqpts, Q.get_weights()) + transform = cell.get_entity_transform(sd-1, entity) + points = transform(Q.get_points()) + weights = numpy.multiply(f_at_qpts, Q.get_weights()) - pt_dict = OrderedDict() - - qpts = mappedQ.pts - qwts = mappedQ.wts - - for k in range(len(qpts)): - pt_cur = tuple(qpts[k]) - pt_dict[pt_cur] = [(qwts[k] * f_at_qpts[k, i, j], (i, j)) - for (i, j) in index_iterator(shp)] + pt_dict = {tuple(pt): [(wt * s1s2T[idx], idx) for idx in index_iterator(shp)] + for pt, wt in zip(points, weights)} super().__init__(cell, shp, pt_dict, {}, nm) @@ -433,15 +424,13 @@ def __init__(self, ref_el, Q, f_at_qpts): sd = ref_el.get_spatial_dimension() - qpts, qwts = Q.get_points(), Q.get_weights() - dpts = qpts - self.dpts = dpts - - dpt_dict = OrderedDict() + points = Q.get_points() + self.dpts = points + weights = numpy.multiply(f_at_qpts, Q.get_weights()) - alphas = [tuple([1 if j == i else 0 for j in range(sd)]) for i in range(sd)] - for j, pt in enumerate(dpts): - dpt_dict[tuple(pt)] = [(qwts[j]*f_at_qpts[j], alphas[i], (i,)) for i in range(sd)] + alphas = tuple(map(tuple, numpy.eye(sd, dtype=int))) + dpt_dict = {tuple(pt): [(wt, alphas[i], (i,)) for i in range(sd)] + for pt, wt in zip(points, weights)} super().__init__(ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfDivergence") @@ -453,24 +442,21 @@ class IntegralMomentOfTensorDivergence(Functional): def __init__(self, ref_el, Q, f_at_qpts): self.f_at_qpts = f_at_qpts self.Q = Q - qpts, qwts = Q.get_points(), Q.get_weights() - nqp = len(qpts) - dpts = qpts - self.dpts = dpts + points = Q.get_points() + self.dpts = points sd = ref_el.get_spatial_dimension() + shp = (sd, sd) assert len(f_at_qpts.shape) == 2 assert f_at_qpts.shape[0] == sd - assert f_at_qpts.shape[1] == nqp + assert f_at_qpts.shape[1] == len(points) + weights = numpy.multiply(f_at_qpts, Q.get_weights()).T - dpt_dict = OrderedDict() + alphas = tuple(map(tuple, numpy.eye(sd, dtype=int))) + dpt_dict = {tuple(pt): [(wt[i], alphas[j], (i, j)) for i, j in index_iterator(shp)] + for pt, wt in zip(points, weights)} - alphas = [tuple([1 if j == i else 0 for j in range(sd)]) for i in range(sd)] - for q, pt in enumerate(dpts): - dpt_dict[tuple(pt)] = [(qwts[q]*f_at_qpts[i, q], alphas[j], (i, j)) for i in range(sd) for j in range(sd)] - - super().__init__(ref_el, tuple(), {}, dpt_dict, - "IntegralMomentOfDivergence") + super().__init__(ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfDivergence") class FrobeniusIntegralMoment(IntegralMoment): @@ -499,11 +485,8 @@ class PointNormalEvaluation(Functional): def __init__(self, ref_el, facet_no, pt): n = ref_el.compute_normal(facet_no) self.n = n - sd = ref_el.get_spatial_dimension() - - pt_dict = {pt: [(n[i], (i,)) for i in range(sd)]} - - shp = (sd,) + shp = n.shape + pt_dict = {pt: [(n[i], (i,)) for i in range(shp[0])]} super().__init__(ref_el, shp, pt_dict, {}, "PointNormalEval") @@ -514,9 +497,8 @@ class PointEdgeTangentEvaluation(Functional): def __init__(self, ref_el, edge_no, pt): t = ref_el.compute_edge_tangent(edge_no) self.t = t - sd = ref_el.get_spatial_dimension() - pt_dict = {pt: [(t[i], (i,)) for i in range(sd)]} - shp = (sd,) + shp = t.shape + pt_dict = {pt: [(t[i], (i,)) for i in range(shp[0])]} super().__init__(ref_el, shp, pt_dict, {}, "PointEdgeTangent") def tostr(self): @@ -539,11 +521,10 @@ def __init__(self, ref_el, Q, P_at_qpts, edge): t = ref_el.compute_edge_tangent(edge) sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(1, edge) - pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) - weights = Q.get_weights() - pt_dict = OrderedDict() - for pt, wgt, phi in zip(pts, weights, P_at_qpts): - pt_dict[pt] = [(wgt*phi*t[i], (i, )) for i in range(sd)] + points = transform(Q.get_points()) + weights = numpy.multiply(P_at_qpts, Q.get_weights()) + pt_dict = {tuple(pt): [(wt*t[i], (i,)) for i in range(sd)] + for pt, wt in zip(points, weights)} super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfEdgeTangentEvaluation") @@ -583,9 +564,9 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): n = ref_el.compute_scaled_normal(facet) sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(sd-1, facet) - pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) + pts = tuple(map(tuple, transform(Q.get_points()))) weights = Q.get_weights() - pt_dict = OrderedDict() + pt_dict = {} for pt, wgt, phi in zip(pts, weights, P_at_qpts): phixn = [phi[1]*n[2] - phi[2]*n[1], phi[2]*n[0] - phi[0]*n[2], @@ -611,12 +592,11 @@ class MonkIntegralMoment(Functional): def __init__(self, ref_el, Q, P_at_qpts, facet): sd = ref_el.get_spatial_dimension() - weights = Q.get_weights() - pt_dict = OrderedDict() transform = ref_el.get_entity_transform(sd-1, facet) - pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) - for pt, wgt, phi in zip(pts, weights, P_at_qpts): - pt_dict[pt] = [(wgt*phi[i], (i, )) for i in range(sd)] + pts = transform(Q.get_points()) + weights = Q.get_weights() * P_at_qpts + pt_dict = {tuple(pt): [(wt[i], (i, )) for i in range(sd)] + for pt, wt in zip(pts, weights)} super().__init__(ref_el, (sd, ), pt_dict, {}, "MonkIntegralMoment") @@ -653,11 +633,10 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): n = ref_el.compute_scaled_normal(facet) sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(sd - 1, facet) - pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) - weights = Q.get_weights() - pt_dict = OrderedDict() - for pt, wgt, phi in zip(pts, weights, P_at_qpts): - pt_dict[pt] = [(wgt*phi*n[i], (i, )) for i in range(sd)] + pts = transform(Q.get_points()) + weights = Q.get_weights() * P_at_qpts + pt_dict = {tuple(pt): [(wt*n[i], (i, )) for i in range(sd)] + for pt, wt in zip(pts, weights)} super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledNormalEvaluation") @@ -672,15 +651,12 @@ class PointwiseInnerProductEvaluation(Functional): correct weights. """ - def __init__(self, ref_el, v, w, p): - sd = ref_el.get_spatial_dimension() - + def __init__(self, ref_el, v, w, pt): wvT = numpy.outer(w, v) + shp = wvT.shape - pt_dict = {p: [(wvT[i][j], (i, j)) - for i, j in index_iterator((sd, sd))]} + pt_dict = {pt: [(wvT[idx], idx) for idx in index_iterator(shp)]} - shp = (sd, sd) super().__init__(ref_el, shp, pt_dict, {}, "PointwiseInnerProductEval") @@ -696,20 +672,15 @@ class TensorBidirectionalMomentInnerProductEvaluation(Functional): """ def __init__(self, ref_el, v, w, Q, f_at_qpts, comp_deg): - sd = ref_el.get_spatial_dimension() - wvT = numpy.outer(w, v) + shp = wvT.shp - qpts, qwts = Q.get_points(), Q.get_weights() + points = Q.get_points() + weights = numpy.multiply(f_at_qpts, Q.get_weights()) - pt_dict = {} - for k, pt in enumerate(map(tuple(qpts))): - pt_dict[pt] = [] - for i, j in index_iterator((sd, sd)): - pt_dict[pt].append((qwts[k] * wvT[i][j] * f_at_qpts[i, j, k]), - (i, j)) + pt_dict = {tuple(pt): [(wt * wvT[idx], idx) for idx in index_iterator(shp)] + for pt, wt in zip(points, weights)} - shp = (sd, sd) super().__init__(ref_el, shp, pt_dict, {}, "TensorBidirectionalMomentInnerProductEvaluation") @@ -728,12 +699,11 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): n = ref_el.compute_scaled_normal(facet) sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(sd - 1, facet) - pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) - weights = Q.get_weights() - pt_dict = OrderedDict() - for pt, wgt, phi in zip(pts, weights, P_at_qpts): - pt_dict[pt] = [(wgt*phi*n[i], (i, )) for i in range(sd)] - super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledNormalEvaluation") + pts = transform(Q.get_points()) + weights = numpy.multiply(P_at_qpts, Q.get_weights()) + pt_dict = {tuple(pt): [(wt*n[i], (i, )) for i in range(sd)] + for pt, wt in zip(pts, weights)} + super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfNormalEvaluation") class IntegralMomentOfTangentialEvaluation(Functional): @@ -752,11 +722,10 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): assert sd == 2 t = ref_el.compute_edge_tangent(facet) transform = ref_el.get_entity_transform(sd - 1, facet) - pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) - weights = Q.get_weights() - pt_dict = OrderedDict() - for pt, wgt, phi in zip(pts, weights, P_at_qpts): - pt_dict[pt] = [(wgt*phi*t[i], (i, )) for i in range(sd)] + points = transform(Q.get_points()) + weights = numpy.multiply(P_at_qpts, Q.get_weights()) + pt_dict = {tuple(pt): [(wt*t[i], (i, )) for i in range(sd)] + for pt, wt in zip(points, weights)} super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledTangentialEvaluation") @@ -773,11 +742,12 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): # scaling on the normal is ok because edge length then weights # the reference element quadrature appropriately n = ref_el.compute_scaled_normal(facet) + nnT = numpy.outer(n, n)/numpy.linalg.norm(n) + shp = nnT.shape sd = ref_el.get_spatial_dimension() transform = ref_el.get_entity_transform(sd - 1, facet) - pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) - weights = Q.get_weights() - pt_dict = OrderedDict() - for pt, wgt, phi in zip(pts, weights, P_at_qpts): - pt_dict[pt] = [(wgt*phi*n[i], (i, )) for i in range(sd)] - super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledNormalEvaluation") + points = transform(Q.get_points()) + weights = numpy.multiply(P_at_qpts, Q.get_weights()) + pt_dict = {tuple(pt): [(wt*nnT[idx], idx) for idx in index_iterator(shp)] + for pt, wt in zip(points, weights)} + super().__init__(ref_el, shp, pt_dict, {}, "IntegralMomentOfNormalNormalEvaluation") diff --git a/FIAT/gauss_legendre.py b/FIAT/gauss_legendre.py index 7962918f6..ff2d9ec82 100644 --- a/FIAT/gauss_legendre.py +++ b/FIAT/gauss_legendre.py @@ -14,4 +14,4 @@ class GaussLegendre(discontinuous_lagrange.DiscontinuousLagrange): """Simplicial discontinuous element with nodes at the (recursive) Gauss-Legendre points.""" def __init__(self, ref_el, degree): - super(GaussLegendre, self).__init__(ref_el, degree, variant="gl") + super().__init__(ref_el, degree, variant="gl") diff --git a/FIAT/gauss_lobatto_legendre.py b/FIAT/gauss_lobatto_legendre.py index e97c7d886..4973a690c 100644 --- a/FIAT/gauss_lobatto_legendre.py +++ b/FIAT/gauss_lobatto_legendre.py @@ -14,4 +14,4 @@ class GaussLobattoLegendre(lagrange.Lagrange): """Simplicial continuous element with nodes at the (recursive) Gauss-Lobatto-Legendre points.""" def __init__(self, ref_el, degree): - super(GaussLobattoLegendre, self).__init__(ref_el, degree, variant="gll", sort_entities=True) + super().__init__(ref_el, degree, variant="gll", sort_entities=True) diff --git a/FIAT/gauss_radau.py b/FIAT/gauss_radau.py index 89191756d..784369b2b 100644 --- a/FIAT/gauss_radau.py +++ b/FIAT/gauss_radau.py @@ -22,7 +22,7 @@ def __init__(self, ref_el, degree, right=True): lr = quadrature.RadauQuadratureLineRule(ref_el, degree+1, right) nodes = [functional.PointEvaluation(ref_el, x) for x in lr.pts] - super(GaussRadauDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class GaussRadau(finite_element.CiarletElement): @@ -33,4 +33,4 @@ def __init__(self, ref_el, degree): poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) dual = GaussRadauDualSet(ref_el, degree) formdegree = ref_el.get_spatial_dimension() # n-form - super(GaussRadau, self).__init__(poly_set, dual, degree, formdegree) + super().__init__(poly_set, dual, degree, formdegree) diff --git a/FIAT/guzman_neilan.py b/FIAT/guzman_neilan.py new file mode 100644 index 000000000..ec72be5a6 --- /dev/null +++ b/FIAT/guzman_neilan.py @@ -0,0 +1,147 @@ +# This file is part of FIAT (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later +# +# Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024 + +# This is not quite Guzman-Neilan, but it has 2*dim*(dim+1) dofs and includes +# dim**2-1 extra constraint functionals. The first (dim+1)**2 basis functions +# are the reference element bfs, but the extra dim**2-1 are used in the +# transformation theory. + +from FIAT import finite_element, polynomial_set, expansions +from FIAT.bernardi_raugel import ExtendedBernardiRaugelSpace, BernardiRaugelDualSet, BernardiRaugel +from FIAT.brezzi_douglas_marini import BrezziDouglasMarini +from FIAT.macro import AlfeldSplit +from FIAT.quadrature_schemes import create_quadrature +from itertools import chain + +import numpy +import math + + +def ExtendedGuzmanNeilanSpace(ref_el, subdegree, reduced=False): + """Return a basis for the extended Guzman-Neilan space.""" + sd = ref_el.get_spatial_dimension() + ref_complex = AlfeldSplit(ref_el) + C0 = polynomial_set.ONPolynomialSet(ref_complex, sd, shape=(sd,), scale=1, variant="bubble") + B = take_interior_bubbles(C0) + if sd > 2: + B = modified_bubble_subspace(B) + + if reduced: + BR = BernardiRaugel(ref_el, sd, subdegree=subdegree).get_nodal_basis() + reduced_dim = BR.get_num_members() - (sd-1) * (sd+1) + BR = BR.take(list(range(reduced_dim))) + else: + BR = ExtendedBernardiRaugelSpace(ref_el, subdegree) + GN = constant_div_projection(BR, C0, B) + return GN + + +class GuzmanNeilan(finite_element.CiarletElement): + """The Guzman-Neilan extended element.""" + def __init__(self, ref_el, degree=None, subdegree=1): + sd = ref_el.get_spatial_dimension() + if degree is None: + degree = sd + if degree != sd: + raise ValueError("Guzman-Neilan only defined for degree = dim") + poly_set = ExtendedGuzmanNeilanSpace(ref_el, subdegree) + dual = BernardiRaugelDualSet(ref_el, degree, subdegree=subdegree) + formdegree = sd - 1 # (n-1)-form + super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") + + +def inner(v, u, qwts): + """Compute the L2 inner product from tabulation arrays and quadrature weights""" + return numpy.tensordot(numpy.multiply(v, qwts), u, + axes=(range(1, v.ndim), range(1, u.ndim))) + + +def div(U): + """Compute the divergence from tabulation dict.""" + return sum(U[k][:, k.index(1), :] for k in U if sum(k) == 1) + + +def take_interior_bubbles(C0): + """Take the interior bubbles from a vector-valued C0 PolynomialSet.""" + ref_complex = C0.get_reference_element() + sd = ref_complex.get_spatial_dimension() + dimC0 = C0.get_num_members() // sd + entity_ids = expansions.polynomial_entity_ids(ref_complex, C0.degree, continuity="C0") + ids = [i + j * dimC0 + for dim in range(sd+1) + for f in sorted(ref_complex.get_interior_facets(dim)) + for i in entity_ids[dim][f] + for j in range(sd)] + return C0.take(ids) + + +def modified_bubble_subspace(B): + """Construct the interior bubble space M_k(K^r) from (Guzman and Neilan, 2019).""" + ref_complex = B.get_reference_element() + sd = ref_complex.get_spatial_dimension() + degree = B.degree + rule = create_quadrature(ref_complex, 2*degree) + qpts, qwts = rule.get_points(), rule.get_weights() + + # tabulate the linear hat function associated with the barycenter + hat = B.take([0]) + hat_at_qpts = hat.tabulate(qpts)[(0,)*sd][0, 0] + + # tabulate the BDM facet functions + ref_el = ref_complex.get_parent() + BDM = BrezziDouglasMarini(ref_el, degree-1) + entity_dofs = BDM.entity_dofs() + facet_dofs = list(range(BDM.space_dimension() - len(entity_dofs[sd][0]))) + BDM_facet = BDM.get_nodal_basis().take(facet_dofs) + phis = BDM_facet.tabulate(qpts)[(0,)*sd] + + # tabulate the bubbles = hat ** (degree - k) * BDMk_facet + bubbles = [numpy.eye(sd)[:, :, None] * hat_at_qpts[None, None, :] ** degree] + for k in range(1, degree): + dimPk = math.comb(k + sd-1, sd-1) + idsPk = list(chain.from_iterable(entity_dofs[sd-1][f][:dimPk] + for f in entity_dofs[sd-1])) + bubbles.append(numpy.multiply(phis[idsPk], hat_at_qpts ** (degree-k))) + bubbles = numpy.concatenate(bubbles, axis=0) + + # store the bubbles into a PolynomialSet via L2 projection + v = B.tabulate(qpts)[(0,) * sd] + coeffs = numpy.linalg.solve(inner(v, v, qwts), inner(v, bubbles, qwts)) + coeffs = numpy.tensordot(coeffs, B.get_coeffs(), axes=(0, 0)) + M = polynomial_set.PolynomialSet(ref_complex, degree, degree, + B.get_expansion_set(), coeffs) + return M + + +def constant_div_projection(BR, C0, M): + """Project the BR space into the space of C0 polynomials with constant divergence.""" + ref_complex = C0.get_reference_element() + sd = ref_complex.get_spatial_dimension() + degree = C0.degree + rule = create_quadrature(ref_complex, 2*degree) + qpts, qwts = rule.get_points(), rule.get_weights() + + # Take the test space for the divergence in L2 \ R + Q = polynomial_set.ONPolynomialSet(ref_complex, degree-1) + Q = Q.take(list(range(1, Q.get_num_members()))) + P = Q.tabulate(qpts)[(0,)*sd] + P -= numpy.dot(P, qwts)[:, None] / sum(qwts) + + U = M.tabulate(qpts, 1) + X = BR.tabulate(qpts, 1) + # Invert the divergence + B = inner(P, div(U), qwts) + g = inner(P, div(X), qwts) + w = numpy.linalg.solve(B, g) + + # Add correction to BR bubbles + v = C0.tabulate(qpts)[(0,)*sd] + coeffs = numpy.linalg.solve(inner(v, v, qwts), inner(v, X[(0,)*sd], qwts)) + coeffs = coeffs.T.reshape(BR.get_num_members(), sd, -1) + coeffs -= numpy.tensordot(w, M.get_coeffs(), axes=(0, 0)) + GN = polynomial_set.PolynomialSet(ref_complex, degree, degree, + C0.get_expansion_set(), coeffs) + return GN diff --git a/FIAT/hct.py b/FIAT/hct.py index 7fe0c146a..b09741cdc 100644 --- a/FIAT/hct.py +++ b/FIAT/hct.py @@ -41,20 +41,19 @@ def __init__(self, ref_complex, degree, reduced=False): entity_ids[0][v].extend(range(cur, len(nodes))) k = 2 if reduced else degree - 3 - rline = ufc_simplex(1) - Q = create_quadrature(rline, degree-1+k) + facet = ufc_simplex(1) + Q = create_quadrature(facet, degree-1+k) qpts = Q.get_points() + xref = 2.0 * qpts - 1.0 if reduced: - x, = qpts.T - f_at_qpts = eval_jacobi(0, 0, k, 2.0*x - 1) + f_at_qpts = eval_jacobi(0, 0, k, xref[:, 0]) for e in sorted(top[1]): cur = len(nodes) nodes.append(IntegralMomentOfNormalDerivative(ref_el, e, Q, f_at_qpts)) entity_ids[1][e].extend(range(cur, len(nodes))) else: - x = 2.0*qpts - 1 - phis = eval_jacobi_batch(1, 1, k, x) - dphis = eval_jacobi_deriv_batch(1, 1, k, x) + phis = eval_jacobi_batch(1, 1, k, xref) + dphis = eval_jacobi_deriv_batch(1, 1, k, xref) for e in sorted(top[1]): Q_mapped = FacetQuadratureRule(ref_el, 1, e, Q) scale = 2 / Q_mapped.jacobian_determinant() @@ -73,7 +72,7 @@ def __init__(self, ref_complex, degree, reduced=False): nodes.extend(IntegralMoment(ref_el, Q, phi * scale) for phi in phis) entity_ids[sd][0] = list(range(cur, len(nodes))) - super(HCTDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class HsiehCloughTocher(finite_element.CiarletElement): @@ -85,4 +84,4 @@ def __init__(self, ref_el, degree=3, reduced=False): ref_complex = macro.AlfeldSplit(ref_el) dual = HCTDualSet(ref_complex, degree, reduced=reduced) poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1, vorder=degree-1, variant="bubble") - super(HsiehCloughTocher, self).__init__(poly_set, dual, degree) + super().__init__(poly_set, dual, degree) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index 0f9ee66c6..a68509f79 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -25,7 +25,7 @@ class TraceError(Exception): or the gradient of a trace element.""" def __init__(self, msg): - super(TraceError, self).__init__(msg) + super().__init__(msg) self.msg = msg @@ -114,9 +114,9 @@ def __init__(self, ref_el, degree): # Degree of the element deg = max([e.degree() for e in dg_elements.values()]) - super(HDivTrace, self).__init__(ref_el, dual, order=deg, - formdegree=facet_sd, - mapping="affine") + super().__init__(ref_el, dual, order=deg, + formdegree=facet_sd, + mapping="affine") # Set up facet elements self.dg_elements = dg_elements diff --git a/FIAT/hellan_herrmann_johnson.py b/FIAT/hellan_herrmann_johnson.py index d8153152f..8218b00f6 100644 --- a/FIAT/hellan_herrmann_johnson.py +++ b/FIAT/hellan_herrmann_johnson.py @@ -39,7 +39,7 @@ def __init__(self, cell, degree): dofs.extend(_dofs) dof_ids[dim] = _dof_ids - super(HellanHerrmannJohnsonDual, self).__init__(dofs, cell, dof_ids) + super().__init__(dofs, cell, dof_ids) @staticmethod def _generate_edge_dofs(cell, degree, offset): @@ -93,5 +93,4 @@ def __init__(self, cell, degree): # mapping under affine transformation mapping = "double contravariant piola" - super(HellanHerrmannJohnson, self).__init__(Ps, Ls, degree, - mapping=mapping) + super().__init__(Ps, Ls, degree, mapping=mapping) diff --git a/FIAT/hermite.py b/FIAT/hermite.py index e1368fea8..79c83ff97 100644 --- a/FIAT/hermite.py +++ b/FIAT/hermite.py @@ -64,7 +64,7 @@ def __init__(self, ref_el): for facet in top[dim]: entity_ids[dim][facet] = [] - super(CubicHermiteDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class CubicHermite(finite_element.CiarletElement): @@ -75,4 +75,4 @@ def __init__(self, ref_el, deg=3): poly_set = polynomial_set.ONPolynomialSet(ref_el, 3) dual = CubicHermiteDualSet(ref_el) - super(CubicHermite, self).__init__(poly_set, dual, 3) + super().__init__(poly_set, dual, 3) diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py index 4120106e8..fcabaf536 100644 --- a/FIAT/hierarchical.py +++ b/FIAT/hierarchical.py @@ -49,7 +49,7 @@ def __init__(self, ref_el, degree, codim=0): entity_ids[dim][entity] = list(range(cur, len(nodes))) entity_permutations[dim][entity] = perms - super(LegendreDual, self).__init__(nodes, ref_el, entity_ids, entity_permutations) + super().__init__(nodes, ref_el, entity_ids, entity_permutations) class Legendre(finite_element.CiarletElement): @@ -60,7 +60,7 @@ def __new__(cls, ref_el, degree, variant=None): if splitting is None: # FIXME P0 on the split requires implementing SplitSimplicialComplex.symmetry_group_size() return P0.P0(ref_el) - return super(Legendre, cls).__new__(cls) + return super().__new__(cls) def __init__(self, ref_el, degree, variant=None): splitting, _ = parse_lagrange_variant(variant, integral=True) @@ -69,7 +69,7 @@ def __init__(self, ref_el, degree, variant=None): poly_set = ONPolynomialSet(ref_el, degree) dual = LegendreDual(ref_el, degree) formdegree = ref_el.get_spatial_dimension() # n-form - super(Legendre, self).__init__(poly_set, dual, degree, formdegree) + super().__init__(poly_set, dual, degree, formdegree) class IntegratedLegendreDual(dual_set.DualSet): @@ -107,7 +107,7 @@ def __init__(self, ref_el, degree): entity_ids[dim][entity] = list(range(cur, len(nodes))) entity_permutations[dim][entity] = perms - super(IntegratedLegendreDual, self).__init__(nodes, ref_el, entity_ids, entity_permutations) + super().__init__(nodes, ref_el, entity_ids, entity_permutations) def make_reference_duals(self, ref_el, degree): Q = create_quadrature(ref_el, 2 * degree) @@ -140,4 +140,4 @@ def __init__(self, ref_el, degree, variant=None): poly_set = ONPolynomialSet(ref_el, degree, variant="bubble") dual = IntegratedLegendreDual(ref_el, degree) formdegree = 0 # 0-form - super(IntegratedLegendre, self).__init__(poly_set, dual, degree, formdegree) + super().__init__(poly_set, dual, degree, formdegree) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index 7f3d9c3d8..de4683bf4 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -18,6 +18,7 @@ def __init__(self, ref_complex, degree, variant=None): nodes = [] # Face dofs: bidirectional (nn and nt) Legendre moments + R = numpy.array([[0, 1], [-1, 0]]) dim = sd - 1 ref_facet = ref_el.construct_subelement(dim) Qref = create_quadrature(ref_facet, 2*degree) @@ -26,12 +27,11 @@ def __init__(self, ref_complex, degree, variant=None): for facet in sorted(top[dim]): cur = len(nodes) Q = FacetQuadratureRule(ref_el, dim, facet, Qref) - Jdet = Q.jacobian_determinant() - tangents = ref_el.compute_tangents(dim, facet) - normal = ref_el.compute_normal(facet) - normal /= numpy.linalg.norm(normal) - scaled_normal = normal * Jdet - uvecs = (scaled_normal, *tangents) + thats = ref_el.compute_tangents(dim, facet) + nhat = numpy.dot(R, *thats) if sd == 2 else numpy.cross(*thats) + normal = nhat / Q.jacobian_determinant() + + uvecs = (nhat, *thats) comps = [numpy.outer(normal, uvec) for uvec in uvecs] nodes.extend(FrobeniusIntegralMoment(ref_el, Q, comp[:, :, None] * phi[None, None, :]) for phi in phis for comp in comps) @@ -48,7 +48,7 @@ def __init__(self, ref_complex, degree, variant=None): entity_ids[sd][0].extend(range(cur, len(nodes))) - super(JohnsonMercierDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class JohnsonMercier(finite_element.CiarletElement): @@ -59,5 +59,4 @@ def __init__(self, ref_el, degree=1, variant=None): poly_set = macro.HDivSymPolynomialSet(ref_complex, degree) dual = JohnsonMercierDualSet(ref_complex, degree, variant=variant) mapping = "double contravariant piola" - super(JohnsonMercier, self).__init__(poly_set, dual, degree, - mapping=mapping) + super().__init__(poly_set, dual, degree, mapping=mapping) diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index 51a996154..eead0cfa6 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -50,7 +50,7 @@ def __init__(self, ref_el, degree, point_variant="equispaced", sort_entities=Fal pts_cur = ref_el.make_points(dim, entity, degree, variant=point_variant) nodes.extend(functional.PointEvaluation(ref_el, x) for x in pts_cur) entity_ids[dim][entity] = list(range(cur, len(nodes))) - super(LagrangeDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) + super().__init__(nodes, ref_el, entity_ids, entity_permutations) class Lagrange(finite_element.CiarletElement): @@ -86,4 +86,4 @@ def __init__(self, ref_el, degree, variant="equispaced", sort_entities=False): poly_variant = "bubble" if ref_el.is_macrocell() else None poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, variant=poly_variant) formdegree = 0 # 0-form - super(Lagrange, self).__init__(poly_set, dual, degree, formdegree) + super().__init__(poly_set, dual, degree, formdegree) diff --git a/FIAT/macro.py b/FIAT/macro.py index 486eda57b..f48298c70 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -1,4 +1,3 @@ -import copy from itertools import chain, combinations import numpy @@ -74,8 +73,9 @@ def make_topology(sd, num_verts, edges): facet = topology[dim][entity] facet_verts = set(facet) for v in range(min(facet)): - if facet_verts < adjacency[v]: + if (facet_verts < adjacency[v]): entities.append((v, *facet)) + topology[dim+1] = dict(enumerate(sorted(entities))) return topology @@ -88,7 +88,10 @@ class SplitSimplicialComplex(SimplicialComplex): :arg topology: The topology of the simplicial complex. """ def __init__(self, parent, vertices, topology): - self._parent = parent + self._parent_complex = parent + while parent.get_parent(): + parent = parent.get_parent() + self._parent_simplex = parent bary = xy_to_bary(parent.get_vertices(), vertices) parent_top = parent.get_topology() @@ -147,7 +150,7 @@ def __init__(self, parent, vertices, topology): for dim in sorted(child_to_parent)} self._interior_facets = interior_facets - super(SplitSimplicialComplex, self).__init__(parent.shape, vertices, topology) + super().__init__(parent.shape, vertices, topology) def get_child_to_parent(self): """Maps split complex facet tuple to its parent entity tuple.""" @@ -180,53 +183,16 @@ def construct_subelement(self, dimension): :arg dimension: subentity dimension (integer) """ - return self._parent.construct_subelement(dimension) + return self.get_parent().construct_subelement(dimension) def is_macrocell(self): return True def get_parent(self): - return self._parent - - -class AlfeldSplit(SplitSimplicialComplex): - """Splits a simplicial complex by connecting subcell vertices to their - barycenter. - """ - def __init__(self, ref_el): - sd = ref_el.get_spatial_dimension() - top = ref_el.get_topology() - # Keep old facets, respecting the old numbering - new_topology = copy.deepcopy(top) - # Discard the cell interiors - new_topology[sd] = {} - new_verts = list(ref_el.get_vertices()) - - for cell in top[sd]: - # Append the barycenter as the new vertex - new_verts.extend(ref_el.make_points(sd, cell, sd+1)) - new_vert_id = len(new_topology[0]) - new_topology[0][new_vert_id] = (new_vert_id,) + return self._parent_simplex - # Append new facets by adding the barycenter to old facets - for dim in range(1, sd + 1): - cur = len(new_topology[dim]) - for entity, ids in top[dim-1].items(): - if set(ids) < set(top[sd][cell]): - new_topology[dim][cur] = ids + (new_vert_id,) - cur = cur + 1 - - parent = ref_el.get_parent() or ref_el - super(AlfeldSplit, self).__init__(parent, tuple(new_verts), new_topology) - - def construct_subcomplex(self, dimension): - """Constructs the reference subcomplex of the parent cell subentity - specified by subcomplex dimension. - """ - if dimension == self.get_dimension(): - return self - # Alfeld on facets is just the parent subcomplex - return self._parent.construct_subcomplex(dimension) + def get_parent_complex(self): + return self._parent_complex class IsoSplit(SplitSimplicialComplex): @@ -264,11 +230,10 @@ def __init__(self, ref_el, degree=2, variant=None): edges.append(tuple(sorted((v0, v1)))) new_topology = make_topology(sd, len(new_verts), edges) - parent = ref_el.get_parent() or ref_el - super(IsoSplit, self).__init__(parent, tuple(new_verts), new_topology) + super().__init__(ref_el, tuple(new_verts), new_topology) def construct_subcomplex(self, dimension): - """Constructs the reference subcomplex of the parent cell subentity + """Constructs the reference subcomplex of the parent complex specified by subcomplex dimension. """ if dimension == self.get_dimension(): @@ -282,43 +247,75 @@ def construct_subcomplex(self, dimension): class PowellSabinSplit(SplitSimplicialComplex): - """Splits a simplicial complex by connecting barycenters of subentities to - the barycenter of the superentity. + """Splits a simplicial complex by connecting barycenters of subentities + to the barycenter of the superentities of the given dimension or higher. """ - def __init__(self, ref_el): + def __init__(self, ref_el, dimension=1): + self.split_dimension = dimension sd = ref_el.get_spatial_dimension() top = ref_el.get_topology() - inv_top = invert_cell_topology(top) + connectivity = ref_el.get_connectivity() new_verts = list(ref_el.get_vertices()) - edges = [] - offset = {0: 0} - for dim in range(1, sd+1): - offset[dim] = len(new_verts) + dim = dimension - 1 + simplices = {dim: {entity: [top[dim][entity]] for entity in top[dim]}} + + for dim in range(dimension, sd+1): + simplices[dim] = {} for entity in top[dim]: - bary_id = entity + offset[dim] + bary_id = len(new_verts) new_verts.extend(ref_el.make_points(dim, entity, dim+1)) - # Connect subentity barycenter to the entity barycenter - for subdim in range(dim): - edges.extend((inv_top[subdim][subverts] + offset[subdim], bary_id) - for subverts in combinations(top[dim][entity], subdim+1)) - - new_topology = make_topology(sd, len(new_verts), edges) - parent = ref_el.get_parent() or ref_el - super(PowellSabinSplit, self).__init__(parent, tuple(new_verts), new_topology) + # Connect entity barycenter to every subsimplex on the entity + simplices[dim][entity] = [(*s, bary_id) + for child in connectivity[(dim, dim-1)][entity] + for s in simplices[dim-1][child]] + + simplices = list(chain.from_iterable(simplices[sd].values())) + new_topology = {} + new_topology[0] = {i: (i,) for i in range(len(new_verts))} + for dim in range(1, sd): + facets = chain.from_iterable((combinations(s, dim+1) for s in simplices)) + if dim < self.split_dimension: + # Preserve the numbering of the unsplit entities + facets = chain(top[dim].values(), facets) + unique_facets = dict.fromkeys(facets) + new_topology[dim] = dict(enumerate(unique_facets)) + new_topology[sd] = dict(enumerate(simplices)) + + parent = ref_el if dimension == sd else PowellSabinSplit(ref_el, dimension=dimension+1) + super().__init__(parent, tuple(new_verts), new_topology) def construct_subcomplex(self, dimension): - """Constructs the reference subcomplex of the parent cell subentity + """Constructs the reference subcomplex of the parent complex specified by subcomplex dimension. """ if dimension == self.get_dimension(): return self - ref_el = self.construct_subelement(dimension) - if dimension == 0: - return ref_el + parent = self.get_parent_complex() + subcomplex = parent.construct_subcomplex(dimension) + if dimension < self.split_dimension: + return subcomplex else: # Powell-Sabin on facets is Powell-Sabin - return PowellSabinSplit(ref_el) + return PowellSabinSplit(subcomplex, dimension=self.split_dimension) + + +class AlfeldSplit(PowellSabinSplit): + """Splits a simplicial complex by connecting cell vertices to their + barycenter. + """ + def __init__(self, ref_el): + sd = ref_el.get_spatial_dimension() + super().__init__(ref_el, dimension=sd) + + +class WorseyFarinSplit(PowellSabinSplit): + """Splits a simplicial complex by connecting cell and facet vertices to their + barycenter. This reduces to Powell-Sabin on the triangle, and Alfeld on the interval. + """ + def __init__(self, ref_el): + sd = ref_el.get_spatial_dimension() + super().__init__(ref_el, dimension=sd-1) class PowellSabin12Split(SplitSimplicialComplex): @@ -345,8 +342,9 @@ def __init__(self, ref_el): (3, 4), (3, 5), (3, 6), (3, 7), (3, 8), (3, 9), (4, 7), (4, 8), (5, 7), (5, 9), (6, 8), (6, 9)] - super(PowellSabin12Split, self).__init__( - ref_el, tuple(new_verts), make_topology(2, len(new_verts), edges)) + parent = PowellSabinSplit(ref_el) + new_topology = make_topology(2, len(new_verts), edges) + super().__init__(parent, tuple(new_verts), new_topology) def construct_subcomplex(self, dimension): """Constructs the reference subcomplex of the parent cell subentity @@ -375,8 +373,7 @@ class MacroQuadratureRule(QuadratureRule): def __init__(self, ref_el, Q_ref, parent_facets=None): parent_dim = Q_ref.ref_el.get_spatial_dimension() if parent_facets is not None: - parent_cell = ref_el.get_parent() - parent_to_children = parent_cell.get_parent_to_children() + parent_to_children = ref_el.get_parent_to_children() facets = [] for parent_entity in parent_facets: children = parent_to_children[parent_dim][parent_entity] @@ -393,7 +390,7 @@ def __init__(self, ref_el, Q_ref, parent_facets=None): wts.extend(Q_cur.wts) pts = tuple(pts) wts = tuple(wts) - super(MacroQuadratureRule, self).__init__(ref_el, pts, wts) + super().__init__(ref_el, pts, wts) class CkPolynomialSet(polynomial_set.PolynomialSet): @@ -452,12 +449,13 @@ def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs): num_sv = len([s for s in sig if abs(s) > tol]) coeffs = numpy.dot(vt[num_sv:], coeffs) - if shape != tuple(): + if shape != (): m, n = coeffs.shape - coeffs = coeffs.reshape((m,) + (1,)*len(shape) + (n,)) - coeffs = numpy.tile(coeffs, (1,) + shape + (1,)) + ncomp = numpy.prod(shape) + coeffs = numpy.kron(coeffs, numpy.eye(ncomp)) + coeffs = coeffs.reshape(m*ncomp, *shape, n) - super(CkPolynomialSet, self).__init__(ref_el, degree, degree, expansion_set, coeffs) + super().__init__(ref_el, degree, degree, expansion_set, coeffs) class HDivSymPolynomialSet(polynomial_set.PolynomialSet): @@ -499,9 +497,9 @@ def __init__(self, ref_el, degree, order=0, **kwargs): rows.append(numpy.tensordot(wn, jump, axes=(ax, ax))) if len(rows) > 0: - dual_mat = numpy.row_stack(rows) + dual_mat = numpy.vstack(rows) _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) num_sv = len([s for s in sig if abs(s) > 1.e-10]) coeffs = numpy.tensordot(vt[num_sv:], coeffs, axes=(1, 0)) - super(HDivSymPolynomialSet, self).__init__(ref_el, degree, degree, expansion_set, coeffs) + super().__init__(ref_el, degree, degree, expansion_set, coeffs) diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index c05133022..c92fa06be 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -171,7 +171,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): for d in range(dim) for phi in Phis) entity_ids[dim][0] = list(range(cur, len(nodes))) - super(NedelecDual, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class Nedelec(finite_element.CiarletElement): @@ -205,5 +205,4 @@ def __init__(self, ref_el, degree, variant=None): raise Exception("Not implemented") dual = NedelecDual(ref_el, degree, variant, interpolant_deg) formdegree = 1 # 1-form - super(Nedelec, self).__init__(poly_set, dual, degree, formdegree, - mapping="covariant piola") + super().__init__(poly_set, dual, degree, formdegree, mapping="covariant piola") diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 6e5269721..481dcf64e 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -52,7 +52,7 @@ def __init__(self, cell, degree, variant, interpolant_deg): # Define degrees of freedom (dofs, ids) = self.generate_degrees_of_freedom(cell, degree, variant, interpolant_deg) # Call init of super-class - super(NedelecSecondKindDual, self).__init__(dofs, cell, ids) + super().__init__(dofs, cell, ids) def generate_degrees_of_freedom(self, cell, degree, variant, interpolant_deg): "Generate dofs and geometry-to-dof maps (ids)." @@ -214,4 +214,4 @@ def __init__(self, cell, degree, variant=None): mapping = "covariant piola" # Call init of super-class - super(NedelecSecondKind, self).__init__(Ps, Ls, degree, formdegree, mapping=mapping) + super().__init__(Ps, Ls, degree, formdegree, mapping=mapping) diff --git a/FIAT/powell_sabin.py b/FIAT/powell_sabin.py index 2b0ea19cb..e9b7c4f3f 100644 --- a/FIAT/powell_sabin.py +++ b/FIAT/powell_sabin.py @@ -37,8 +37,7 @@ def __init__(self, ref_complex, degree=2): nodes.extend(PointDerivative(ref_el, pt, alpha) for alpha in alphas) entity_ids[0][v].extend(range(cur, len(nodes))) - super(QuadraticPowellSabin6DualSet, self).__init__( - nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class QuadraticPowellSabin6(finite_element.CiarletElement): @@ -52,7 +51,7 @@ def __init__(self, ref_el, degree=2): dual = QuadraticPowellSabin6DualSet(ref_complex, degree) poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1) - super(QuadraticPowellSabin6, self).__init__(poly_set, dual, degree) + super().__init__(poly_set, dual, degree) class QuadraticPowellSabin12DualSet(dual_set.DualSet): @@ -90,8 +89,7 @@ def __init__(self, ref_complex, degree=2): nodes.extend(IntegralMomentOfNormalDerivative(ref_el, e, Q, phi) for phi in phis) entity_ids[1][e].extend(range(cur, len(nodes))) - super(QuadraticPowellSabin12DualSet, self).__init__( - nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class QuadraticPowellSabin12(finite_element.CiarletElement): @@ -105,4 +103,4 @@ def __init__(self, ref_el, degree=2): dual = QuadraticPowellSabin12DualSet(ref_complex, degree) poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1) - super(QuadraticPowellSabin12, self).__init__(poly_set, dual, degree) + super().__init__(poly_set, dual, degree) diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py index e8d26aacc..61702c54c 100644 --- a/FIAT/quadrature.py +++ b/FIAT/quadrature.py @@ -21,14 +21,16 @@ def pseudo_determinant(A): def map_quadrature(pts_ref, wts_ref, source_cell, target_cell, jacobian=False): """Map quadrature points and weights defined on source_cell to target_cell. """ + while source_cell.get_parent(): + source_cell = source_cell.get_parent() A, b = reference_element.make_affine_mapping(source_cell.get_vertices(), target_cell.get_vertices()) + if len(pts_ref.shape) != 2: + pts_ref = pts_ref.reshape(-1, A.shape[1]) scale = pseudo_determinant(A) + pts = numpy.dot(pts_ref, A.T) + pts = numpy.add(pts, b, out=pts) wts = scale * wts_ref - if pts_ref.size == 0: - pts = b[None, :] - else: - pts = numpy.dot(pts_ref.reshape((-1, A.shape[1])), A.T) + b[None, :] # return immutable types pts = tuple(map(tuple, pts)) @@ -92,7 +94,7 @@ class GaussLegendreQuadratureLineRule(GaussJacobiQuadratureLineRule): The quadrature rule uses m points for a degree of precision of 2m-1. """ def __init__(self, ref_el, m): - super(GaussLegendreQuadratureLineRule, self).__init__(ref_el, m) + super().__init__(ref_el, m) class RadauQuadratureLineRule(QuadratureRule): diff --git a/FIAT/quadrature_element.py b/FIAT/quadrature_element.py index d71e59c7c..d61eff908 100644 --- a/FIAT/quadrature_element.py +++ b/FIAT/quadrature_element.py @@ -32,7 +32,7 @@ def __init__(self, ref_el, points, weights=None): # Construct the dual set dual = DualSet(nodes, ref_el, entity_dofs) - super(QuadratureElement, self).__init__(ref_el, dual, order=None) + super().__init__(ref_el, dual, order=None) self._points = points # save the quadrature points & weights self._weights = weights diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index f33ce937f..db0058d12 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -15,7 +15,7 @@ def RTSpace(ref_el, degree): - """Constructs a basis for the the Raviart-Thomas space + """Constructs a basis for the Raviart-Thomas space (P_{degree-1})^d + P_{degree-1} x""" sd = ref_el.get_spatial_dimension() @@ -116,7 +116,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): for pt in pts) entity_ids[sd][0] = list(range(cur, len(nodes))) - super(RTDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) class RaviartThomas(finite_element.CiarletElement): @@ -146,5 +146,4 @@ def __init__(self, ref_el, degree, variant=None): poly_set = RTSpace(ref_el, degree) dual = RTDualSet(ref_el, degree, variant, interpolant_deg) formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form - super(RaviartThomas, self).__init__(poly_set, dual, degree, formdegree, - mapping="contravariant piola") + super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 1e750bbed..f09751578 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -177,12 +177,6 @@ def _key(self): # Default: only type matters return None - def __eq__(self, other): - return type(self) is type(other) and self._key() == other._key() - - def __ne__(self, other): - return not self.__eq__(other) - def __hash__(self): return hash((type(self), self._key())) @@ -268,21 +262,56 @@ def is_simplex(self): def is_macrocell(self): return False + def get_interior_facets(self, dim): + """Return the interior facets this cell is a split and () otherwise.""" + return () + def get_parent(self): """Return the parent cell if this cell is a split and None otherwise.""" return None + def get_parent_complex(self): + """Return the parent complex if this cell is a split and None otherwise.""" + return None + + def is_parent(self, other, strict=False): + """Return whether this cell is the parent of the other cell.""" + parent = other + if strict: + parent = parent.get_parent_complex() + while parent is not None: + if self == parent: + return True + parent = parent.get_parent_complex() + return False + + def __eq__(self, other): + if self is other: + return True + A, B = self.get_vertices(), other.get_vertices() + if not (len(A) == len(B) and numpy.allclose(A, B)): + return False + atop = self.get_topology() + btop = other.get_topology() + for dim in atop: + if set(atop[dim].values()) != set(btop[dim].values()): + return False + return True + + def __ne__(self, other): + return not self.__eq__(other) + def __gt__(self, other): - return self.get_parent() == other + return other.is_parent(self, strict=True) def __lt__(self, other): - return self == other.get_parent() + return self.is_parent(other, strict=True) def __ge__(self, other): - return self > other or self == other + return other.is_parent(self, strict=False) def __le__(self, other): - return self < other or self == other + return self.is_parent(other, strict=False) class SimplicialComplex(Cell): @@ -297,7 +326,7 @@ def __init__(self, shape, vertices, topology): for entity in topology[dim]: assert len(topology[dim][entity]) == dim + 1 - super(SimplicialComplex, self).__init__(shape, vertices, topology) + super().__init__(shape, vertices, topology) def compute_normal(self, facet_i, cell=None): """Returns the unit normal vector to facet i of codimension 1.""" @@ -386,7 +415,7 @@ def compute_normalized_tangents(self, dim, i): of dimension dim. Returns a (possibly empty) list. These tangents are normalized to have unit length.""" ts = self.compute_tangents(dim, i) - ts /= numpy.linalg.norm(ts, axis=0)[None, :] + ts /= numpy.linalg.norm(ts, axis=1)[:, None] return ts def compute_edge_tangent(self, edge_i): @@ -539,6 +568,11 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): out = numpy.dot(points, A.T) return numpy.add(out, b, out=out) + def compute_bubble(self, points, entity=None): + """Returns the lowest-order bubble on an entity evaluated at the given + points on the entity.""" + return numpy.prod(self.compute_barycentric_coordinates(points, entity), axis=1) + def distance_to_point_l1(self, points, entity=None, rescale=False): # noqa: D301 """Get the L1 distance (aka 'manhatten', 'taxicab' or rectilinear @@ -784,7 +818,7 @@ class Point(Simplex): def __init__(self): verts = ((),) topology = {0: {0: (0,)}} - super(Point, self).__init__(POINT, verts, topology) + super().__init__(POINT, verts, topology) def construct_subelement(self, dimension): """Constructs the reference element of a cell subentity @@ -804,7 +838,7 @@ def __init__(self): edges = {0: (0, 1)} topology = {0: {0: (0,), 1: (1,)}, 1: edges} - super(DefaultLine, self).__init__(LINE, verts, topology) + super().__init__(LINE, verts, topology) class UFCInterval(UFCSimplex): @@ -815,7 +849,7 @@ def __init__(self): edges = {0: (0, 1)} topology = {0: {0: (0,), 1: (1,)}, 1: edges} - super(UFCInterval, self).__init__(LINE, verts, topology) + super().__init__(LINE, verts, topology) class DefaultTriangle(DefaultSimplex): @@ -830,7 +864,7 @@ def __init__(self): faces = {0: (0, 1, 2)} topology = {0: {0: (0,), 1: (1,), 2: (2,)}, 1: edges, 2: faces} - super(DefaultTriangle, self).__init__(TRIANGLE, verts, topology) + super().__init__(TRIANGLE, verts, topology) class UFCTriangle(UFCSimplex): @@ -843,7 +877,7 @@ def __init__(self): faces = {0: (0, 1, 2)} topology = {0: {0: (0,), 1: (1,), 2: (2,)}, 1: edges, 2: faces} - super(UFCTriangle, self).__init__(TRIANGLE, verts, topology) + super().__init__(TRIANGLE, verts, topology) def compute_normal(self, i): "UFC consistent normal" @@ -863,7 +897,7 @@ def __init__(self): faces = {0: (0, 1, 2)} topology = {0: {0: (0,), 1: (1,), 2: (2,)}, 1: edges, 2: faces} - super(IntrepidTriangle, self).__init__(TRIANGLE, verts, topology) + super().__init__(TRIANGLE, verts, topology) def get_facet_element(self): # I think the UFC interval is equivalent to what the @@ -894,7 +928,7 @@ def __init__(self): 3: (0, 1, 2)} tets = {0: (0, 1, 2, 3)} topology = {0: vs, 1: edges, 2: faces, 3: tets} - super(DefaultTetrahedron, self).__init__(TETRAHEDRON, verts, topology) + super().__init__(TETRAHEDRON, verts, topology) class IntrepidTetrahedron(Simplex): @@ -919,7 +953,7 @@ def __init__(self): 3: (0, 2, 1)} tets = {0: (0, 1, 2, 3)} topology = {0: vs, 1: edges, 2: faces, 3: tets} - super(IntrepidTetrahedron, self).__init__(TETRAHEDRON, verts, topology) + super().__init__(TETRAHEDRON, verts, topology) def get_facet_element(self): return IntrepidTriangle() @@ -947,7 +981,7 @@ def __init__(self): 3: (0, 1, 2)} tets = {0: (0, 1, 2, 3)} topology = {0: vs, 1: edges, 2: faces, 3: tets} - super(UFCTetrahedron, self).__init__(TETRAHEDRON, verts, topology) + super().__init__(TETRAHEDRON, verts, topology) def compute_normal(self, i): "UFC consistent normals." @@ -982,7 +1016,7 @@ def __init__(self, *cells): topology[dim] = dict(enumerate(topology[dim][key] for key in sorted(topology[dim]))) - super(TensorProductCell, self).__init__(TENSORPRODUCT, vertices, topology) + super().__init__(TENSORPRODUCT, vertices, topology) self.cells = tuple(cells) def _key(self): @@ -1181,7 +1215,7 @@ def __init__(self): verts = product.get_vertices() topology = flatten_entities(pt) - super(UFCQuadrilateral, self).__init__(QUADRILATERAL, verts, topology) + super().__init__(QUADRILATERAL, verts, topology) self.product = product self.unflattening_map = compute_unflattening_map(pt) @@ -1284,7 +1318,7 @@ def __init__(self): verts = product.get_vertices() topology = flatten_entities(pt) - super(UFCHexahedron, self).__init__(HEXAHEDRON, verts, topology) + super().__init__(HEXAHEDRON, verts, topology) self.product = product self.unflattening_map = compute_unflattening_map(pt) diff --git a/FIAT/regge.py b/FIAT/regge.py index e6f90e6f6..e9292f663 100644 --- a/FIAT/regge.py +++ b/FIAT/regge.py @@ -42,7 +42,7 @@ def __init__(self, cell, degree): dofs.extend(_dofs) dof_ids[dim] = _dof_ids - super(ReggeDual, self).__init__(dofs, cell, dof_ids) + super().__init__(dofs, cell, dof_ids) @staticmethod def _generate_dofs(cell, entity_dim, degree, offset): @@ -92,4 +92,4 @@ def __init__(self, cell, degree): # mapping under affine transformation mapping = "double covariant piola" - super(Regge, self).__init__(Ps, Ls, degree, mapping=mapping) + super().__init__(Ps, Ls, degree, mapping=mapping) diff --git a/FIAT/restricted.py b/FIAT/restricted.py index 26c82d913..45d73c099 100644 --- a/FIAT/restricted.py +++ b/FIAT/restricted.py @@ -12,24 +12,19 @@ class RestrictedDualSet(DualSet): """Restrict the given DualSet to the specified list of dofs.""" def __init__(self, dual, indices): + indices = list(sorted(indices)) ref_el = dual.get_reference_element() nodes_old = dual.get_nodes() - dof_counter = 0 entity_ids = {} nodes = [] for d, entities in dual.get_entity_ids().items(): entity_ids[d] = {} for entity, dofs in entities.items(): - entity_ids[d][entity] = [] - for dof in dofs: - if dof not in indices: - continue - entity_ids[d][entity].append(dof_counter) - dof_counter += 1 - nodes.append(nodes_old[dof]) - assert dof_counter == len(indices) + entity_ids[d][entity] = [indices.index(dof) + for dof in dofs if dof in indices] + nodes = [nodes_old[i] for i in indices] self._dual = dual - super(RestrictedDualSet, self).__init__(nodes, ref_el, entity_ids) + super().__init__(nodes, ref_el, entity_ids) def get_indices(self, restriction_domain, take_closure=True): """Return the list of dofs with support on a given restriction domain. @@ -74,4 +69,4 @@ def __init__(self, element, indices=None, restriction_domain=None, take_closure= assert all(e_mapping == mapping_new[0] for e_mapping in mapping_new) # Call constructor of CiarletElement - super(RestrictedElement, self).__init__(poly_set, dual, 0, element.get_formdegree(), mapping_new[0]) + super().__init__(poly_set, dual, 0, element.get_formdegree(), mapping_new[0]) diff --git a/FIAT/serendipity.py b/FIAT/serendipity.py index 4f783dd58..94b7930db 100644 --- a/FIAT/serendipity.py +++ b/FIAT/serendipity.py @@ -121,7 +121,7 @@ def __init__(self, ref_el, degree): assert len(s_list) == cur formdegree = 0 - super(Serendipity, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree) + super().__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree) self.basis = {(0,)*dim: Array(s_list)} polynomials, extra_vars = _replace_numbers_with_symbols(Array(s_list)) diff --git a/FIAT/tensor_product.py b/FIAT/tensor_product.py index 90c59c91d..9d2671d2e 100644 --- a/FIAT/tensor_product.py +++ b/FIAT/tensor_product.py @@ -206,7 +206,7 @@ def __init__(self, A, B): dual = dual_set.DualSet(nodes, ref_el, entity_ids) - super(TensorProductElement, self).__init__(ref_el, dual, order, formdegree, mapping) + super().__init__(ref_el, dual, order, formdegree, mapping) # Set up constituent elements self.A = A self.B = B @@ -383,7 +383,7 @@ def __init__(self, element): flat_entity_ids = flatten_entities(entity_ids) dual = DualSet(nodes, ref_el, flat_entity_ids) - super(FlattenedDimensions, self).__init__(ref_el, dual, element.get_order(), element.get_formdegree(), element._mapping) + super().__init__(ref_el, dual, element.get_order(), element.get_formdegree(), element._mapping) self.element = element # Construct unflattening map for passing correct values to tabulate() diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index bed985f5d..5e2cf264d 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -42,9 +42,16 @@ from FIAT.tensor_product import TensorProductElement # noqa: F401 from FIAT.tensor_product import FlattenedDimensions # noqa: F401 from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401 +from FIAT.bernardi_raugel import BernardiRaugel # noqa: F401 from FIAT.argyris import Argyris # noqa: F401 from FIAT.hermite import CubicHermite # noqa: F401 from FIAT.morley import Morley # noqa: F401 +from FIAT.hct import HsiehCloughTocher # noqa: F401 +from FIAT.alfeld_sorokina import AlfeldSorokina # noqa: F401 +from FIAT.arnold_qin import ArnoldQin # noqa: F401 +from FIAT.christiansen_hu import ChristiansenHu # noqa: F401 +from FIAT.guzman_neilan import GuzmanNeilan # noqa: F401 +from FIAT.johnson_mercier import JohnsonMercier # noqa: F401 from FIAT.bubble import Bubble from FIAT.enriched import EnrichedElement # noqa: F401 from FIAT.nodal_enriched import NodalEnrichedElement @@ -298,6 +305,25 @@ def __init__(self, a, b): "CubicHermite(T)", "CubicHermite(S)", "Morley(T)", + "BernardiRaugel(T)", + "BernardiRaugel(S)", + + # Macroelements + "Lagrange(T, 1, 'iso')", + "Lagrange(T, 1, 'alfeld')", + "Lagrange(T, 2, 'alfeld')", + "DiscontinuousLagrange(T, 1, 'alfeld')", + "HsiehCloughTocher(T)", + "JohnsonMercier(T)", + "JohnsonMercier(S)", + "AlfeldSorokina(T)", + "AlfeldSorokina(S)", + "ArnoldQin(T, reduced=False)", + "ArnoldQin(T, reduced=True)", + "ChristiansenHu(T)", + "ChristiansenHu(S)", + "GuzmanNeilan(T)", + "GuzmanNeilan(S)", # MixedElement made of nodal elements should be nodal, but its API # is currently just broken. @@ -351,7 +377,7 @@ def test_nodality(element): # Fetch primal and dual basis poly_set = element.get_nodal_basis() dual_set = element.get_dual_set() - assert poly_set.get_reference_element() == dual_set.get_reference_element() + assert poly_set.get_reference_element() >= dual_set.get_reference_element() # Get coeffs of primal and dual bases w.r.t. expansion set coeffs_poly = poly_set.get_coeffs() diff --git a/test/unit/test_hct.py b/test/unit/test_hct.py index 27b67560d..4dbb96024 100644 --- a/test/unit/test_hct.py +++ b/test/unit/test_hct.py @@ -8,7 +8,9 @@ @pytest.fixture def cell(): - return ufc_simplex(2) + K = ufc_simplex(2) + K.vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84)) + return K @pytest.mark.parametrize("reduced", (False, True)) diff --git a/test/unit/test_macro.py b/test/unit/test_macro.py index 4a6646716..f1aa683c3 100644 --- a/test/unit/test_macro.py +++ b/test/unit/test_macro.py @@ -2,7 +2,7 @@ import numpy import pytest from FIAT import DiscontinuousLagrange, Lagrange, Legendre, P0 -from FIAT.macro import AlfeldSplit, IsoSplit, CkPolynomialSet +from FIAT.macro import AlfeldSplit, IsoSplit, PowellSabinSplit, CkPolynomialSet from FIAT.quadrature_schemes import create_quadrature from FIAT.reference_element import ufc_simplex from FIAT.expansions import polynomial_entity_ids, polynomial_cell_node_map @@ -10,9 +10,9 @@ from FIAT.barycentric_interpolation import get_lagrange_points -@pytest.fixture(params=("I", "T", "S")) +@pytest.fixture(params=(1, 2, 3), ids=("I", "T", "S")) def cell(request): - dim = {"I": 1, "T": 2, "S": 3}[request.param] + dim = request.param return ufc_simplex(dim) @@ -129,6 +129,21 @@ def test_macro_lagrange(variant, degree, split, cell): assert numpy.allclose(fe.V, V) +def test_powell_sabin(cell): + dim = cell.get_spatial_dimension() + A = AlfeldSplit(cell) + assert A > cell + + PS = PowellSabinSplit(cell, dim) + assert PS == A + + for split_dim in range(1, dim): + PS = PowellSabinSplit(cell, split_dim) + assert PS > A + assert PS > cell + assert len(PS.get_topology()[dim]) == math.factorial(dim+1) // math.factorial(split_dim) + + def make_mass_matrix(fe, order=0): sd = fe.ref_el.get_spatial_dimension() Q = create_quadrature(fe.ref_complex, 2*fe.degree()) diff --git a/test/unit/test_stokes_complex.py b/test/unit/test_stokes_complex.py new file mode 100644 index 000000000..1c20b0e75 --- /dev/null +++ b/test/unit/test_stokes_complex.py @@ -0,0 +1,144 @@ +import pytest +import numpy + +from FIAT import HsiehCloughTocher as HCT +from FIAT import AlfeldSorokina as AS +from FIAT import ArnoldQin as AQ +from FIAT import Lagrange as CG +from FIAT import DiscontinuousLagrange as DG +from FIAT.restricted import RestrictedElement +from FIAT.reference_element import ufc_simplex + +from FIAT.macro import CkPolynomialSet +from FIAT.alfeld_sorokina import AlfeldSorokinaSpace +from FIAT.arnold_qin import ArnoldQinSpace +from FIAT.christiansen_hu import ChristiansenHuSpace +from FIAT.guzman_neilan import ExtendedGuzmanNeilanSpace + + +T = ufc_simplex(2) +S = ufc_simplex(3) + + +@pytest.mark.parametrize("cell", (T, S)) +@pytest.mark.parametrize("degree", (2, 3, 4)) +def test_AlfeldSorokinaSpace(cell, degree): + # Test that the divergence of the Alfeld-Sorokina space is spanned by a C0 basis + V = AlfeldSorokinaSpace(cell, degree) + A = V.get_reference_element() + top = A.get_topology() + sd = A.get_spatial_dimension() + + pts = [] + for dim in top: + for entity in top[dim]: + pts.extend(A.make_points(dim, entity, degree)) + + V_tab = V.tabulate(pts, 1) + V_div = sum(V_tab[alpha][:, alpha.index(1), :] + for alpha in V_tab if sum(alpha) == 1) + + C0 = CkPolynomialSet(A, degree-1, order=0, variant="bubble") + C0_tab = C0.tabulate(pts)[(0,)*sd] + _, residual, *_ = numpy.linalg.lstsq(C0_tab.T, V_div.T) + assert numpy.allclose(residual, 0) + + +@pytest.mark.parametrize("cell", (T,)) +@pytest.mark.parametrize("reduced", (False, True)) +def test_hct_stokes_complex(cell, reduced): + # Test that we have the lowest-order discrete Stokes complex, verifying + # that the range of the exterior derivative of each space is contained by + # the next space in the sequence + H2 = HCT(T, reduced=reduced) + A = H2.get_reference_complex() + if reduced: + # HCT-red(3) -curl-> AQ-red(2) -div-> DG(0) + H2 = RestrictedElement(H2, restriction_domain="vertex") + H1 = RestrictedElement(AQ(T, reduced=reduced), indices=list(range(9))) + L2 = DG(cell, 0) + else: + # HCT(3) -curl-> AS(2) -div-> CG(1, Alfeld) + H1 = AS(cell) + L2 = CG(A, 1) + + pts = [] + top = A.get_topology() + for dim in top: + for entity in top[dim]: + pts.extend(A.make_points(dim, entity, 4)) + + H2tab = H2.tabulate(1, pts) + H1tab = H1.tabulate(1, pts) + L2tab = L2.tabulate(0, pts) + + L2val = L2tab[(0, 0)] + H1val = H1tab[(0, 0)] + H1div = sum(H1tab[alpha][:, alpha.index(1), :] for alpha in H1tab if sum(alpha) == 1) + H2curl = numpy.stack([H2tab[(0, 1)], -H2tab[(1, 0)]], axis=1) + + H2dim = H2.space_dimension() + H1dim = H1.space_dimension() + _, residual, *_ = numpy.linalg.lstsq(H1val.reshape(H1dim, -1).T, H2curl.reshape(H2dim, -1).T) + assert numpy.allclose(residual, 0) + + _, residual, *_ = numpy.linalg.lstsq(L2val.T, H1div.T) + assert numpy.allclose(residual, 0) + + +@pytest.mark.parametrize("cell", (T, S)) +@pytest.mark.parametrize("family", ("AQ", "CH", "GN")) +def test_minimal_stokes_space(cell, family): + # Test that the C0 Stokes space is spanned by a C0 basis + # Also test that its divergence is constant + sd = cell.get_spatial_dimension() + if family == "GN": + degree = 1 + space = ExtendedGuzmanNeilanSpace + elif family == "CH": + degree = 1 + space = ChristiansenHuSpace + elif family == "AQ": + if sd != 2: + return + degree = 2 + space = ArnoldQinSpace + + W = space(cell, degree) + V = space(cell, degree, reduced=True) + Wdim = W.get_num_members() + Vdim = V.get_num_members() + K = W.get_reference_element() + sd = K.get_spatial_dimension() + top = K.get_topology() + + pts = [] + for dim in top: + for entity in top[dim]: + pts.extend(K.make_points(dim, entity, degree)) + + C0 = CkPolynomialSet(K, sd, order=0, variant="bubble") + C0_tab = C0.tabulate(pts) + Wtab = W.tabulate(pts, 1) + Vtab = V.tabulate(pts, 1) + z = (0,)*sd + for tab in (Vtab, Wtab): + # Test that the space is full rank + _, sig, _ = numpy.linalg.svd(tab[z].reshape(-1, sd*len(pts)).T, full_matrices=True) + assert all(sig > 1E-10) + + # Test that the space is C0 + for k in range(sd): + _, residual, *_ = numpy.linalg.lstsq(C0_tab[z].T, tab[z][:, k, :].T) + assert numpy.allclose(residual, 0) + + # Test that divergence is in P0 + div = sum(tab[alpha][:, alpha.index(1), :] + for alpha in tab if sum(alpha) == 1)[:Vdim] + assert numpy.allclose(div, div[:, 0][:, None]) + + # Test that the full space includes the reduced space + assert Wdim > Vdim + _, residual, *_ = numpy.linalg.lstsq(Wtab[z].reshape(Wdim, -1).T, + Vtab[z].reshape(Vdim, -1).T) + assert numpy.allclose(residual, 0) From bf9130f2f689968dada5697dbacb6a5ddd9a3f8e Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 14 Oct 2024 09:07:23 +0100 Subject: [PATCH 04/13] Fix intersphinx mapping --- doc/sphinx/source/conf.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/doc/sphinx/source/conf.py b/doc/sphinx/source/conf.py index cf9642bc2..f7c7b2797 100644 --- a/doc/sphinx/source/conf.py +++ b/doc/sphinx/source/conf.py @@ -263,7 +263,10 @@ # Example configuration for intersphinx: refer to the Python standard library. -intersphinx_mapping = {'http://docs.python.org/': None} +intersphinx_mapping = { + 'numpy': ('https://numpy.org/doc/stable/', None), + 'python': ('https://docs.python.org/3/', None), +} def run_apidoc(_): From 30afc82f384360d6c83ebb932866e65cc3d422a3 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 14 Oct 2024 09:12:21 +0100 Subject: [PATCH 05/13] add recursivenodes --- doc/sphinx/source/conf.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/doc/sphinx/source/conf.py b/doc/sphinx/source/conf.py index f7c7b2797..f97faaf30 100644 --- a/doc/sphinx/source/conf.py +++ b/doc/sphinx/source/conf.py @@ -262,8 +262,9 @@ #texinfo_no_detailmenu = False -# Example configuration for intersphinx: refer to the Python standard library. +# Configuration for intersphinx intersphinx_mapping = { + 'recursivenodes': ('https://tisaac.gitlab.io/recursivenodes/', None), 'numpy': ('https://numpy.org/doc/stable/', None), 'python': ('https://docs.python.org/3/', None), } From 11be0ebb631993924aabf62b5ea6fc0f86a48ad4 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 16 Oct 2024 16:48:39 +0100 Subject: [PATCH 06/13] NodalEnriched macroelements (#86) * NodalEnriched macroelements --- FIAT/__init__.py | 6 +- FIAT/alfeld_sorokina.py | 48 ++++-------- FIAT/arnold_qin.py | 12 +-- FIAT/bernardi_raugel.py | 110 ++++++++++++++------------- FIAT/christiansen_hu.py | 4 +- FIAT/dual_set.py | 34 ++++++--- FIAT/expansions.py | 5 ++ FIAT/finite_element.py | 2 +- FIAT/guzman_neilan.py | 126 ++++++++++++++++++++++++------- FIAT/hierarchical.py | 34 ++++----- FIAT/nodal_enriched.py | 73 ++++++++++-------- FIAT/polynomial_set.py | 10 +-- test/unit/test_fiat.py | 12 ++- test/unit/test_stokes_complex.py | 102 ++++++++++++++++++------- 14 files changed, 355 insertions(+), 223 deletions(-) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 65afac645..117606fd7 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -13,7 +13,7 @@ from FIAT.hct import HsiehCloughTocher from FIAT.alfeld_sorokina import AlfeldSorokina from FIAT.arnold_qin import ArnoldQin -from FIAT.guzman_neilan import GuzmanNeilan +from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanH1div from FIAT.christiansen_hu import ChristiansenHu from FIAT.johnson_mercier import JohnsonMercier from FIAT.brezzi_douglas_marini import BrezziDouglasMarini @@ -88,7 +88,9 @@ "Alfeld-Sorokina": AlfeldSorokina, "Arnold-Qin": ArnoldQin, "Christiansen-Hu": ChristiansenHu, - "Guzman-Neilan": GuzmanNeilan, + "Guzman-Neilan 1st kind H1": GuzmanNeilanFirstKindH1, + "Guzman-Neilan 2nd kind H1": GuzmanNeilanSecondKindH1, + "Guzman-Neilan H1(div)": GuzmanNeilanH1div, "Johnson-Mercier": JohnsonMercier, "Lagrange": Lagrange, "Kong-Mulder-Veldhuizen": KongMulderVeldhuizen, diff --git a/FIAT/alfeld_sorokina.py b/FIAT/alfeld_sorokina.py index 13ef25bce..5d3244b78 100644 --- a/FIAT/alfeld_sorokina.py +++ b/FIAT/alfeld_sorokina.py @@ -7,18 +7,16 @@ # Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024 from FIAT import finite_element, dual_set, polynomial_set -from FIAT.functional import ComponentPointEvaluation, PointDivergence, IntegralMoment +from FIAT.functional import ComponentPointEvaluation, PointDivergence from FIAT.quadrature_schemes import create_quadrature -from FIAT.quadrature import FacetQuadratureRule from FIAT.macro import CkPolynomialSet, AlfeldSplit import numpy def AlfeldSorokinaSpace(ref_el, degree): - """Return a vector-valued C^0 PolynomialSet on an Alfeld split whose - divergence is also C^0. This works on any simplex and for all - polynomial degrees.""" + """Return a vector-valued C0 PolynomialSet on an Alfeld split with C0 + divergence. This works on any simplex and for all polynomial degrees.""" ref_complex = AlfeldSplit(ref_el) sd = ref_complex.get_spatial_dimension() C0 = CkPolynomialSet(ref_complex, degree, order=0, shape=(sd,), variant="bubble") @@ -51,53 +49,35 @@ def AlfeldSorokinaSpace(ref_el, degree): class AlfeldSorokinaDualSet(dual_set.DualSet): - def __init__(self, ref_el, degree, variant=None): + def __init__(self, ref_el, degree): if degree != 2: - raise NotImplementedError("Alfeld-Sorokina only defined for degree = 2") - if variant is None: - variant = "integral" - if variant not in {"integral", "point"}: - raise ValueError(f"Invalid variant {variant}") + raise NotImplementedError(f"{type(self).__name__} only defined for degree = 2") top = ref_el.get_topology() sd = ref_el.get_spatial_dimension() entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} nodes = [] - dims = (0, 1) if variant == "point" else (0,) - for dim in dims: + for dim in sorted(top): for entity in sorted(top[dim]): cur = len(nodes) + + dpts = ref_el.make_points(dim, entity, degree-1) + nodes.extend(PointDivergence(ref_el, pt) for pt in dpts) + pts = ref_el.make_points(dim, entity, degree) - if dim == 0: - pt, = pts - nodes.append(PointDivergence(ref_el, pt)) nodes.extend(ComponentPointEvaluation(ref_el, k, (sd,), pt) for pt in pts for k in range(sd)) entity_ids[dim][entity].extend(range(cur, len(nodes))) - if variant == "integral": - # Edge moments against quadratic Legendre (mean-free bubble) - dim = 1 - facet = ref_el.construct_subelement(dim) - Q = create_quadrature(facet, degree+dim+1) - f_at_qpts = facet.compute_bubble(Q.get_points()) - f_at_qpts -= numpy.dot(f_at_qpts, Q.get_weights()) / facet.volume() - for entity in sorted(top[dim]): - cur = len(nodes) - Q_mapped = FacetQuadratureRule(ref_el, dim, entity, Q) - detJ = Q_mapped.jacobian_determinant() - phi = f_at_qpts / detJ - nodes.extend(IntegralMoment(ref_el, Q_mapped, phi, comp=k, shp=(sd,)) - for k in range(sd)) - entity_ids[dim][entity].extend(range(cur, len(nodes))) - super().__init__(nodes, ref_el, entity_ids) class AlfeldSorokina(finite_element.CiarletElement): - """The Alfeld-Sorokina C^0 quadratic macroelement with C^0 divergence. This element - belongs to a Stokes complex, and is paired with Lagrange(ref_el, 1, variant="alfeld").""" + """The Alfeld-Sorokina C0 quadratic macroelement with C0 divergence. + + This element belongs to a Stokes complex, and is paired with CG1(Alfeld). + """ def __init__(self, ref_el, degree=2): dual = AlfeldSorokinaDualSet(ref_el, degree) poly_set = AlfeldSorokinaSpace(ref_el, degree) diff --git a/FIAT/arnold_qin.py b/FIAT/arnold_qin.py index c2a421dd5..1c5794da8 100644 --- a/FIAT/arnold_qin.py +++ b/FIAT/arnold_qin.py @@ -16,8 +16,8 @@ def ArnoldQinSpace(ref_el, degree, reduced=False): """Return a basis for the Arnold-Qin space - curl(HCT-red) + P_0 x if reduced = True, and - curl(HCT) + P_0 x if reduced = False.""" + curl(HCT-red) + P0 x if reduced = True, and + curl(HCT) + P0 x if reduced = False.""" if ref_el.get_shape() != TRIANGLE: raise ValueError("Arnold-Qin only defined on triangles") if degree != 2: @@ -56,16 +56,16 @@ def ArnoldQinSpace(ref_el, degree, reduced=False): class ArnoldQin(finite_element.CiarletElement): - """The Arnold-Qin C^0(Alfeld) quadratic macroelement with divergence in P0. + """The Arnold-Qin C0(Alfeld) quadratic macroelement with divergence in P0. This element belongs to a Stokes complex, and is paired with unsplit DG0.""" def __init__(self, ref_el, degree=2, reduced=False): poly_set = ArnoldQinSpace(ref_el, degree) if reduced: - subdegree = 1 + order = 1 mapping = "contravariant piola" else: - subdegree = degree + order = degree mapping = "affine" - dual = BernardiRaugelDualSet(ref_el, degree, subdegree) + dual = BernardiRaugelDualSet(ref_el, order, degree=degree) formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) diff --git a/FIAT/bernardi_raugel.py b/FIAT/bernardi_raugel.py index 68e8720f2..6809dfdd0 100644 --- a/FIAT/bernardi_raugel.py +++ b/FIAT/bernardi_raugel.py @@ -11,92 +11,98 @@ from FIAT import finite_element, dual_set, polynomial_set, expansions from FIAT.functional import ComponentPointEvaluation, FrobeniusIntegralMoment -from FIAT.quadrature_schemes import create_quadrature +from FIAT.hierarchical import make_dual_bubbles from FIAT.quadrature import FacetQuadratureRule + import numpy +import math -def ExtendedBernardiRaugelSpace(ref_el, subdegree): - r"""Return a basis for the extended Bernardi-Raugel space. - P_k^d + (P_{dim} \ P_{dim-1})^d""" +def BernardiRaugelSpace(ref_el, order): + """Return a basis for the extended Bernardi-Raugel space: (Pk + FacetBubble)^d.""" sd = ref_el.get_spatial_dimension() - if subdegree >= sd: - raise ValueError("The Bernardi-Raugel space is only defined for subdegree < dim") + if order > sd: + raise ValueError("The Bernardi-Raugel space is only defined for order <= dim") Pd = polynomial_set.ONPolynomialSet(ref_el, sd, shape=(sd,), scale=1, variant="bubble") dimPd = expansions.polynomial_dimension(ref_el, sd, continuity="C0") entity_ids = expansions.polynomial_entity_ids(ref_el, sd, continuity="C0") + + slices = {dim: slice(math.comb(order-1, dim)) for dim in range(order)} + slices[sd-1] = slice(None) ids = [i + j * dimPd - for dim in (*tuple(range(subdegree)), sd-1) + for dim in slices for f in sorted(entity_ids[dim]) - for i in entity_ids[dim][f] + for i in entity_ids[dim][f][slices[dim]] for j in range(sd)] return Pd.take(ids) class BernardiRaugelDualSet(dual_set.DualSet): """The Bernardi-Raugel dual set.""" - def __init__(self, ref_complex, degree, subdegree=1, reduced=False): - ref_el = ref_complex.get_parent() or ref_complex + def __init__(self, ref_el, order=1, degree=None, reduced=False, ref_complex=None): + if ref_complex is None: + ref_complex = ref_el sd = ref_el.get_spatial_dimension() - if subdegree > sd: - raise ValueError("The Bernardi-Raugel dual is only defined for subdegree <= dim") + if degree is None: + degree = sd + if order > sd: + raise ValueError(f"{type(self).__name__} is only defined for order <= dim") top = ref_el.get_topology() entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} - # Point evaluation at lattice points nodes = [] - for dim in range(subdegree): - for entity in sorted(top[dim]): - cur = len(nodes) - pts = ref_el.make_points(dim, entity, subdegree) - nodes.extend(ComponentPointEvaluation(ref_el, comp, (sd,), pt) - for pt in pts for comp in range(sd)) - entity_ids[dim][entity].extend(range(cur, len(nodes))) - - if subdegree < sd: - # Face moments of normal/tangential components against mean-free bubbles - facet = ref_complex.construct_subcomplex(sd-1) - Q = create_quadrature(facet, 2*degree) - if degree == 1 and facet.is_macrocell(): - P = polynomial_set.ONPolynomialSet(facet, degree, scale=1, variant="bubble") - f_at_qpts = P.tabulate(Q.get_points())[(0,)*(sd-1)][-1] - else: - ref_facet = facet.get_parent() or facet - f_at_qpts = ref_facet.compute_bubble(Q.get_points()) - f_at_qpts -= numpy.dot(f_at_qpts, Q.get_weights()) / facet.volume() - - Qs = {f: FacetQuadratureRule(ref_el, sd-1, f, Q) - for f in sorted(top[sd-1])} - - thats = {f: ref_el.compute_tangents(sd-1, f) - for f in sorted(top[sd-1])} + if order > 0: + # Point evaluation at lattice points + for dim in sorted(top): + for entity in sorted(top[dim]): + cur = len(nodes) + pts = ref_el.make_points(dim, entity, order) + nodes.extend(ComponentPointEvaluation(ref_el, comp, (sd,), pt) + for pt in pts for comp in range(sd)) + entity_ids[dim][entity].extend(range(cur, len(nodes))) + + if order < sd: + # Face moments of normal/tangential components against dual bubbles + ref_facet = ref_complex.construct_subcomplex(sd-1) + codim = sd-1 if degree == 1 and ref_facet.is_macrocell() else 0 + Q, phis = make_dual_bubbles(ref_facet, degree, codim=codim) + f_at_qpts = phis[-1] + if codim != 0: + f_at_qpts -= numpy.dot(f_at_qpts, Q.get_weights()) / ref_facet.volume() + + interior_facets = ref_el.get_interior_facets(sd-1) or () + facets = list(set(top[sd-1]) - set(interior_facets)) + + Qs = {f: FacetQuadratureRule(ref_el, sd-1, f, Q) for f in facets} + thats = {f: ref_el.compute_tangents(sd-1, f) for f in facets} R = numpy.array([[0, 1], [-1, 0]]) ndir = 1 if reduced else sd for i in range(ndir): - for f, Q_mapped in Qs.items(): + for f in sorted(facets): cur = len(nodes) if i == 0: udir = numpy.dot(R, *thats[f]) if sd == 2 else numpy.cross(*thats[f]) else: udir = thats[f][i-1] - detJ = Q_mapped.jacobian_determinant() + detJ = Qs[f].jacobian_determinant() phi_at_qpts = udir[:, None] * f_at_qpts[None, :] / detJ - nodes.append(FrobeniusIntegralMoment(ref_el, Q_mapped, phi_at_qpts)) + nodes.append(FrobeniusIntegralMoment(ref_el, Qs[f], phi_at_qpts)) entity_ids[sd-1][f].extend(range(cur, len(nodes))) - super().__init__(nodes, ref_el, entity_ids) class BernardiRaugel(finite_element.CiarletElement): - """The Bernardi-Raugel extended element.""" - def __init__(self, ref_el, degree=None, subdegree=1): - sd = ref_el.get_spatial_dimension() - if degree is None: - degree = sd - if degree != sd: - raise ValueError("Bernardi-Raugel only defined for degree = dim") - poly_set = ExtendedBernardiRaugelSpace(ref_el, subdegree) - dual = BernardiRaugelDualSet(ref_el, degree, subdegree=subdegree) - formdegree = sd - 1 # (n-1)-form + """The Bernardi-Raugel (extended) element. + + This element does not belong to a Stokes complex, but can be paired with + DG_{k-1}. This pair is inf-sup stable, but only weakly divergence-free. + """ + def __init__(self, ref_el, order=1): + degree = ref_el.get_spatial_dimension() + if order >= degree: + raise ValueError(f"{type(self).__name__} only defined for order < dim") + poly_set = BernardiRaugelSpace(ref_el, order) + dual = BernardiRaugelDualSet(ref_el, order, degree=degree) + formdegree = 0 super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/christiansen_hu.py b/FIAT/christiansen_hu.py index d9d5b33e0..476ea914d 100644 --- a/FIAT/christiansen_hu.py +++ b/FIAT/christiansen_hu.py @@ -40,7 +40,7 @@ def ChristiansenHuSpace(ref_el, degree, reduced=False): if not reduced: # Compute the primal basis via Vandermonde and extract the facet bubbles - dual = BernardiRaugelDualSet(ref_complex, degree, reduced=True) + dual = BernardiRaugelDualSet(ref_el, degree, degree=degree, ref_complex=ref_complex, reduced=True) dualmat = dual.to_riesz(C0) V = numpy.tensordot(dualmat, coeffs, axes=((1, 2), (1, 2))) coeffs = numpy.tensordot(numpy.linalg.inv(V.T), coeffs, axes=(-1, 0)) @@ -73,6 +73,6 @@ def __init__(self, ref_el, degree=1): raise ValueError("Christiansen-Hu only defined for degree = 1") poly_set = ChristiansenHuSpace(ref_el, degree) ref_complex = poly_set.get_reference_element() - dual = BernardiRaugelDualSet(ref_complex, degree) + dual = BernardiRaugelDualSet(ref_el, degree, degree=degree, ref_complex=ref_complex) formdegree = ref_el.get_spatial_dimension() - 1 # (n-1)-form super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index d6d0de44b..438ffb9af 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -253,7 +253,7 @@ def make_entity_closure_ids(ref_el, entity_ids): def lexsort_nodes(ref_el, nodes, entity=None, offset=0): """Sort PointEvaluation nodes in lexicographical ordering.""" - if len(nodes) > 1 and all(isinstance(node, functional.PointEvaluation) for node in nodes): + if len(nodes) > 1: pts = [] for node in nodes: pt, = node.get_point_dict() @@ -270,17 +270,29 @@ 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_nodes = [] parent_ids = {} parent_permutations = None - parent_to_children = ref_el.get_parent_to_children() - 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]: - parent_nodes.extend(nodes[i] for i in entity_ids[child_dim][child_entity]) - ids = lexsort_nodes(parent_cell, parent_nodes[cur:], entity=(dim, entity), offset=cur) - parent_ids[dim][entity] = ids + + if all(isinstance(node, functional.PointEvaluation) for node in nodes): + # 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]: + parent_nodes.extend(nodes[i] for i in entity_ids[child_dim][child_entity]) + ids = lexsort_nodes(parent_cell, parent_nodes[cur:], entity=(dim, entity), offset=cur) + parent_ids[dim][entity] = ids + else: + # 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]) + return parent_nodes, parent_cell, parent_ids, parent_permutations diff --git a/FIAT/expansions.py b/FIAT/expansions.py index ac0eea173..7fa0ad887 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -512,6 +512,11 @@ def tabulate_jet(self, n, pts, order=1): data.append(vr.transpose((r, r+1) + tuple(range(r)))) return data + def __eq__(self, other): + return (type(self) is type(other) and + self.ref_el == other.ref_el and + self.continuity == other.continuity) + class PointExpansionSet(ExpansionSet): """Evaluates the point basis on a point reference element.""" diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index 20d05f142..6f4ff24d5 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -132,7 +132,7 @@ class CiarletElement(FiniteElement): def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref_complex=None): ref_el = dual.get_reference_element() ref_complex = ref_complex or poly_set.get_reference_element() - super(CiarletElement, self).__init__(ref_el, dual, order, formdegree, mapping, ref_complex) + super().__init__(ref_el, dual, order, formdegree, mapping, ref_complex) # build generalized Vandermonde matrix old_coeffs = poly_set.get_coeffs() diff --git a/FIAT/guzman_neilan.py b/FIAT/guzman_neilan.py index ec72be5a6..0c721a308 100644 --- a/FIAT/guzman_neilan.py +++ b/FIAT/guzman_neilan.py @@ -10,18 +10,33 @@ # transformation theory. from FIAT import finite_element, polynomial_set, expansions -from FIAT.bernardi_raugel import ExtendedBernardiRaugelSpace, BernardiRaugelDualSet, BernardiRaugel +from FIAT.bernardi_raugel import BernardiRaugelSpace, BernardiRaugelDualSet, BernardiRaugel +from FIAT.alfeld_sorokina import AlfeldSorokina from FIAT.brezzi_douglas_marini import BrezziDouglasMarini from FIAT.macro import AlfeldSplit from FIAT.quadrature_schemes import create_quadrature +from FIAT.restricted import RestrictedElement +from FIAT.nodal_enriched import NodalEnrichedElement from itertools import chain import numpy import math -def ExtendedGuzmanNeilanSpace(ref_el, subdegree, reduced=False): - """Return a basis for the extended Guzman-Neilan space.""" +def GuzmanNeilanSpace(ref_el, order, kind=1, reduced=False): + r"""Return a basis for the (extended) Guzman-Neilan H1 space. + + Project the extended Bernardi-Raugel space (Pk + FacetBubble)^d + into C0 Pk(Alfeld)^d with P_{k-1} divergence, preserving its trace. + + :arg ref_el: a simplex + :arg order: the maximal polynomial degree + :kwarg kind: kind = 1 gives Pk^d + GN bubbles, + kind = 2 gives C0 Pk(Alfeld)^d + GN bubbles. + :kwarg reduced: Include tangential bubbles if reduced = False. + + :returns: a PolynomialSet basis for the Guzman-Neilan H1 space. + """ sd = ref_el.get_spatial_dimension() ref_complex = AlfeldSplit(ref_el) C0 = polynomial_set.ONPolynomialSet(ref_complex, sd, shape=(sd,), scale=1, variant="bubble") @@ -30,29 +45,79 @@ def ExtendedGuzmanNeilanSpace(ref_el, subdegree, reduced=False): B = modified_bubble_subspace(B) if reduced: - BR = BernardiRaugel(ref_el, sd, subdegree=subdegree).get_nodal_basis() + BR = BernardiRaugel(ref_el, order).get_nodal_basis() reduced_dim = BR.get_num_members() - (sd-1) * (sd+1) BR = BR.take(list(range(reduced_dim))) else: - BR = ExtendedBernardiRaugelSpace(ref_el, subdegree) + BR = BernardiRaugelSpace(ref_el, order) + GN = constant_div_projection(BR, C0, B) + if kind == 2: + Bk = take_interior_bubbles(C0, order) + GN = polynomial_set.polynomial_set_union_normalized(GN, Bk) return GN -class GuzmanNeilan(finite_element.CiarletElement): - """The Guzman-Neilan extended element.""" - def __init__(self, ref_el, degree=None, subdegree=1): +class GuzmanNeilanH1(finite_element.CiarletElement): + """The Guzman-Neilan H1-conforming (extended) macroelement.""" + def __init__(self, ref_el, order=1, kind=1): sd = ref_el.get_spatial_dimension() - if degree is None: - degree = sd - if degree != sd: - raise ValueError("Guzman-Neilan only defined for degree = dim") - poly_set = ExtendedGuzmanNeilanSpace(ref_el, subdegree) - dual = BernardiRaugelDualSet(ref_el, degree, subdegree=subdegree) + if order >= sd: + raise ValueError(f"{type(self).__name__} is only defined for order < dim") + degree = sd + poly_set = GuzmanNeilanSpace(ref_el, order, kind=kind) + ref_complex = poly_set.get_reference_element() if kind == 2 else ref_el + dual = BernardiRaugelDualSet(ref_complex, order, degree=degree) formdegree = sd - 1 # (n-1)-form super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") +class GuzmanNeilanFirstKindH1(GuzmanNeilanH1): + """The Guzman-Neilan H1-conforming (extended) macroelement of the first kind. + + Reference element: a simplex of any dimension. + Function space: Pk^d + normal facet bubbles with div in P0, with 1 <= k < dim. + Degrees of freedom: evaluation at Pk lattice points, and normal moments on faces. + + This element belongs to a Stokes complex, and is paired with unsplit DG_{k-1}. + """ + def __init__(self, ref_el, order=1): + super().__init__(ref_el, order=order, kind=1) + + +class GuzmanNeilanSecondKindH1(GuzmanNeilanH1): + """The Guzman-Neilan H1-conforming (extended) macroelement of the second kind. + + Reference element: a simplex of any dimension. + Function space: C0 Pk^d(Alfeld) + normal facet bubbles with div in P0, with 1 <= k < dim. + Degrees of freedom: evaluation at Pk(Alfeld) lattice points, and normal moments on faces. + + This element belongs to a Stokes complex, and is paired with DG_{k-1}(Alfeld). + """ + def __init__(self, ref_el, order=1): + super().__init__(ref_el, order=order, kind=2) + + +def GuzmanNeilanH1div(ref_el, degree=2, reduced=False): + """The Guzman-Neilan H1(div)-conforming (extended) macroelement. + + Reference element: a simplex of any dimension. + Function space: C0 P2^d(Alfeld) with C0 P1 divergence + normal facet bubbles with div in P0. + Degrees of freedom: evaluation at P2(Alfeld) lattice points, divergence at P1 lattice points, + and normal moments on faces. + + This element belongs to a Stokes complex, and is paired with CG1(Alfeld). + """ + AS = AlfeldSorokina(ref_el, 2) + if reduced: + AS = RestrictedElement(AS, restriction_domain="vertex") + elif ref_el.get_spatial_dimension() <= 2: + # Quadratic bubbles are already included in 2D + return AS + GN = GuzmanNeilanH1(ref_el, order=0) + return NodalEnrichedElement(AS, GN) + + def inner(v, u, qwts): """Compute the L2 inner product from tabulation arrays and quadrature weights""" return numpy.tensordot(numpy.multiply(v, qwts), u, @@ -64,18 +129,26 @@ def div(U): return sum(U[k][:, k.index(1), :] for k in U if sum(k) == 1) -def take_interior_bubbles(C0): - """Take the interior bubbles from a vector-valued C0 PolynomialSet.""" - ref_complex = C0.get_reference_element() - sd = ref_complex.get_spatial_dimension() - dimC0 = C0.get_num_members() // sd - entity_ids = expansions.polynomial_entity_ids(ref_complex, C0.degree, continuity="C0") - ids = [i + j * dimC0 - for dim in range(sd+1) +def take_interior_bubbles(P, degree=None): + """Extract the interior bubbles up to the given degree from a complete PolynomialSet.""" + ref_complex = P.get_reference_element() + ncomp = numpy.prod(P.get_shape()) + dimPk = P.expansion_set.get_num_members(P.degree) + assert ncomp * dimPk == P.get_num_members() + continuity = P.expansion_set.continuity + entity_ids = expansions.polynomial_entity_ids(ref_complex, P.degree, + continuity=continuity) + if degree is None or degree >= P.degree: + slices = {dim: slice(None) for dim in entity_ids} + else: + slices = {dim: slice(math.comb(degree-1, dim)) for dim in entity_ids} + + ids = [i + j * dimPk + for dim in slices for f in sorted(ref_complex.get_interior_facets(dim)) - for i in entity_ids[dim][f] - for j in range(sd)] - return C0.take(ids) + for i in entity_ids[dim][f][slices[dim]] + for j in range(ncomp)] + return P.take(ids) def modified_bubble_subspace(B): @@ -117,7 +190,8 @@ def modified_bubble_subspace(B): def constant_div_projection(BR, C0, M): - """Project the BR space into the space of C0 polynomials with constant divergence.""" + """Project the BR space into C0 Pk(Alfeld)^d with P_{k-1} divergence. + """ ref_complex = C0.get_reference_element() sd = ref_complex.get_spatial_dimension() degree = C0.degree diff --git a/FIAT/hierarchical.py b/FIAT/hierarchical.py index fcabaf536..e8666572b 100644 --- a/FIAT/hierarchical.py +++ b/FIAT/hierarchical.py @@ -7,7 +7,6 @@ # Written by Pablo D. Brubeck (brubeck@protonmail.com), 2022 import numpy -import scipy from FIAT import finite_element, dual_set, functional, P0 from FIAT.reference_element import symmetric_simplex @@ -18,6 +17,18 @@ from FIAT.check_format_variant import parse_lagrange_variant +def make_dual_bubbles(ref_el, degree, codim=0): + """Tabulate the L2-duals of the hierarchical C0 basis.""" + Q = create_quadrature(ref_el, 2 * degree) + B = make_bubbles(ref_el, degree, codim=codim, scale="orthonormal") + P_at_qpts = B.expansion_set.tabulate(degree, Q.get_points()) + + M = numpy.dot(numpy.multiply(P_at_qpts, Q.get_weights()), P_at_qpts.T) + phis = numpy.linalg.solve(M, P_at_qpts) + phis = numpy.dot(B.get_coeffs(), phis) + return Q, phis + + class LegendreDual(dual_set.DualSet): """The dual basis for Legendre elements.""" def __init__(self, ref_el, degree, codim=0): @@ -94,7 +105,7 @@ def __init__(self, ref_el, degree): continue ref_facet = symmetric_simplex(dim) - Q_ref, phis = self.make_reference_duals(ref_facet, degree) + Q_ref, phis = make_dual_bubbles(ref_facet, degree) for entity in sorted(top[dim]): cur = len(nodes) Q_facet = FacetQuadratureRule(ref_el, dim, entity, Q_ref) @@ -109,25 +120,6 @@ def __init__(self, ref_el, degree): super().__init__(nodes, ref_el, entity_ids, entity_permutations) - def make_reference_duals(self, ref_el, degree): - Q = create_quadrature(ref_el, 2 * degree) - qpts, qwts = Q.get_points(), Q.get_weights() - inner = lambda v, u: numpy.dot(numpy.multiply(v, qwts), u.T) - dim = ref_el.get_spatial_dimension() - - B = make_bubbles(ref_el, degree) - B_table = B.expansion_set.tabulate(degree, qpts) - - P = ONPolynomialSet(ref_el, degree) - P_table = P.tabulate(qpts, 0)[(0,) * dim] - - # TODO sparse LU - V = inner(P_table, B_table) - PLU = scipy.linalg.lu_factor(V) - phis = scipy.linalg.lu_solve(PLU, P_table) - phis = numpy.dot(B.get_coeffs(), phis) - return Q, phis - class IntegratedLegendre(finite_element.CiarletElement): """Simplicial continuous element with integrated Legendre polynomials.""" diff --git a/FIAT/nodal_enriched.py b/FIAT/nodal_enriched.py index b98eac831..1b258e6d6 100644 --- a/FIAT/nodal_enriched.py +++ b/FIAT/nodal_enriched.py @@ -5,7 +5,9 @@ # SPDX-License-Identifier: LGPL-3.0-or-later import numpy as np +import math +from FIAT.expansions import polynomial_entity_ids from FIAT.polynomial_set import PolynomialSet from FIAT.dual_set import DualSet from FIAT.finite_element import CiarletElement @@ -16,8 +18,8 @@ class NodalEnrichedElement(CiarletElement): """NodalEnriched element is a direct sum of a sequence of - finite elements. Dual basis is reorthogonalized to the - primal basis for nodality. + finite elements. Primal basis is reorthogonalized to the + dual basis for nodality. The following is equivalent: * the constructor is well-defined, @@ -36,38 +38,34 @@ def __init__(self, *elements): "of NodalEnrichedElement are nodal") # Extract common data - degree = max(e.get_nodal_basis().get_degree() for e in elements) - embedded_degree = max(e.get_nodal_basis().get_embedded_degree() - for e in elements) + embedded_degrees = [e.get_nodal_basis().get_embedded_degree() for e in elements] + embedded_degree = max(embedded_degrees) + degree = max(e.degree() for e in elements) order = max(e.get_order() for e in elements) formdegree = None if any(e.get_formdegree() is None for e in elements) \ else max(e.get_formdegree() for e in elements) - # LagrangeExpansionSet has fixed degree, ensure we grab the embedding one - elem = next(e for e in elements - if e.get_nodal_basis().get_embedded_degree() == embedded_degree) - ref_el = elem.get_reference_element() + # LagrangeExpansionSet has fixed degree, ensure we grab the highest one + elem = elements[embedded_degrees.index(embedded_degree)] + ref_el = elem.get_reference_complex() expansion_set = elem.get_nodal_basis().get_expansion_set() mapping = elem.mapping()[0] value_shape = elem.value_shape() # Sanity check - assert all(e.get_nodal_basis().get_reference_element() == - ref_el for e in elements) - assert all(e_mapping == mapping for e in elements - for e_mapping in e.mapping()) + assert all(e.get_reference_complex() == ref_el for e in elements) + assert all(set(e.mapping()) == {mapping, } for e in elements) assert all(e.value_shape() == value_shape for e in elements) # Merge polynomial sets if isinstance(expansion_set, LagrangeLineExpansionSet): # Obtain coefficients via interpolation points = expansion_set.get_points() - coeffs = [e.tabulate(0, points)[(0,)] for e in elements] + coeffs = np.vstack([e.tabulate(0, points)[(0,)] for e in elements]) else: - assert all(type(e.get_nodal_basis().get_expansion_set()) == - type(expansion_set) for e in elements) + assert all(e.get_nodal_basis().get_expansion_set() == expansion_set + for e in elements) coeffs = [e.get_coeffs() for e in elements] - - coeffs = _merge_coeffs(coeffs) + coeffs = _merge_coeffs(coeffs, ref_el, embedded_degrees, expansion_set.continuity) poly_set = PolynomialSet(ref_el, degree, embedded_degree, @@ -81,15 +79,18 @@ def __init__(self, *elements): # Merge dual bases nodes = [node for e in elements for node in e.dual_basis()] + ref_el = ref_el.get_parent() or ref_el dual_set = DualSet(nodes, ref_el, entity_ids) # CiarletElement constructor adjusts poly_set coefficients s.t. # dual_set is really dual to poly_set - super(NodalEnrichedElement, self).__init__(poly_set, dual_set, order, - formdegree=formdegree, mapping=mapping) + super().__init__(poly_set, dual_set, order, formdegree=formdegree, mapping=mapping) + +def _merge_coeffs(coeffss, ref_el, degrees, continuity): + # Indices of the hierachical expansion set on each facet + entity_ids = polynomial_entity_ids(ref_el, max(degrees), continuity) -def _merge_coeffs(coeffss): # Number of bases members total_dim = sum(c.shape[0] for c in coeffss) @@ -101,14 +102,26 @@ def _merge_coeffs(coeffss): max_expansion_dim = max(c.shape[-1] for c in coeffss) # Compose new coeffs - shape = (total_dim,) + value_shape + (max_expansion_dim,) + shape = (total_dim, *value_shape, max_expansion_dim) new_coeffs = np.zeros(shape, dtype=coeffss[0].dtype) counter = 0 - for c in coeffss: - dim = c.shape[0] - expansion_dim = c.shape[-1] - new_coeffs[counter:counter+dim, ..., :expansion_dim] = c - counter += dim + for c, degree in zip(coeffss, degrees): + ids = [] + if continuity == "C0": + dims = sorted(entity_ids) + else: + dims = (ref_el.get_spatial_dimension(),) + for dim in dims: + if continuity == "C0": + dimPk = math.comb(degree - 1, dim) + else: + dimPk = math.comb(degree + dim, dim) + for entity in sorted(entity_ids[dim]): + ids.extend(entity_ids[dim][entity][:dimPk]) + + num_members = c.shape[0] + new_coeffs[counter:counter+num_members, ..., ids] = c + counter += num_members assert counter == total_dim return new_coeffs @@ -117,10 +130,10 @@ def _merge_entity_ids(entity_ids, offsets): ret = {} for i, ids in enumerate(entity_ids): for dim in ids: - if not ret.get(dim): + if dim not in ret: ret[dim] = {} for entity in ids[dim]: - if not ret[dim].get(entity): + if entity not in ret[dim]: ret[dim][entity] = [] - ret[dim][entity] += (np.array(ids[dim][entity]) + offsets[i]).tolist() + ret[dim][entity].extend(np.array(ids[dim][entity]) + offsets[i]) return ret diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 732d03892..579ec3614 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -136,8 +136,7 @@ def __init__(self, ref_el, degree, shape=tuple(), **kwargs): coeffs[cur_idx] = 1.0 cur_bf += 1 - super(ONPolynomialSet, self).__init__(ref_el, degree, embedded_degree, - expansion_set, coeffs) + super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs) def project(f, U, Q): @@ -224,14 +223,13 @@ def __init__(self, ref_el, degree, size=None, **kwargs): coeffs[cur_bf, j, i, exp_bf] = 1.0 cur_bf += 1 - super(ONSymTensorPolynomialSet, self).__init__(ref_el, degree, embedded_degree, - expansion_set, coeffs) + super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs) -def make_bubbles(ref_el, degree, codim=0, shape=()): +def make_bubbles(ref_el, degree, codim=0, shape=(), scale="L2 piola"): """Construct a polynomial set with codim bubbles up to the given degree. """ - poly_set = ONPolynomialSet(ref_el, degree, shape=shape, scale="L2 piola", variant="bubble") + poly_set = ONPolynomialSet(ref_el, degree, shape=shape, scale=scale, variant="bubble") entity_ids = expansions.polynomial_entity_ids(ref_el, degree, continuity="C0") sd = ref_el.get_spatial_dimension() dim = sd - codim diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 5e2cf264d..729aa310e 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -50,7 +50,8 @@ from FIAT.alfeld_sorokina import AlfeldSorokina # noqa: F401 from FIAT.arnold_qin import ArnoldQin # noqa: F401 from FIAT.christiansen_hu import ChristiansenHu # noqa: F401 -from FIAT.guzman_neilan import GuzmanNeilan # noqa: F401 +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.enriched import EnrichedElement # noqa: F401 @@ -322,8 +323,13 @@ def __init__(self, a, b): "ArnoldQin(T, reduced=True)", "ChristiansenHu(T)", "ChristiansenHu(S)", - "GuzmanNeilan(T)", - "GuzmanNeilan(S)", + "GuzmanNeilanFirstKindH1(T, 1)", + "GuzmanNeilanFirstKindH1(S, 1)", + "GuzmanNeilanFirstKindH1(S, 2)", + "GuzmanNeilanSecondKindH1(T, 1)", + "GuzmanNeilanSecondKindH1(S, 1)", + "GuzmanNeilanSecondKindH1(S, 2)", + "NodalEnrichedElement(GuzmanNeilanFirstKindH1(S, 0), AlfeldSorokina(S))", # MixedElement made of nodal elements should be nodal, but its API # is currently just broken. diff --git a/test/unit/test_stokes_complex.py b/test/unit/test_stokes_complex.py index 1c20b0e75..d163ae60c 100644 --- a/test/unit/test_stokes_complex.py +++ b/test/unit/test_stokes_complex.py @@ -1,47 +1,89 @@ import pytest import numpy -from FIAT import HsiehCloughTocher as HCT -from FIAT import AlfeldSorokina as AS -from FIAT import ArnoldQin as AQ -from FIAT import Lagrange as CG -from FIAT import DiscontinuousLagrange as DG -from FIAT.restricted import RestrictedElement +from FIAT import (HsiehCloughTocher as HCT, + AlfeldSorokina as AS, + ArnoldQin as AQ, + Lagrange as CG, + DiscontinuousLagrange as DG) from FIAT.reference_element import ufc_simplex +from FIAT.polynomial_set import ONPolynomialSet from FIAT.macro import CkPolynomialSet from FIAT.alfeld_sorokina import AlfeldSorokinaSpace from FIAT.arnold_qin import ArnoldQinSpace from FIAT.christiansen_hu import ChristiansenHuSpace -from FIAT.guzman_neilan import ExtendedGuzmanNeilanSpace +from FIAT.guzman_neilan import GuzmanNeilanSpace, GuzmanNeilanH1div +from FIAT.restricted import RestrictedElement T = ufc_simplex(2) S = ufc_simplex(3) -@pytest.mark.parametrize("cell", (T, S)) -@pytest.mark.parametrize("degree", (2, 3, 4)) -def test_AlfeldSorokinaSpace(cell, degree): - # Test that the divergence of the Alfeld-Sorokina space is spanned by a C0 basis - V = AlfeldSorokinaSpace(cell, degree) +def span_greater_equal(A, B): + # span(A) >= span(B) + dimA = A.shape[0] + dimB = B.shape[0] + _, residual, *_ = numpy.linalg.lstsq(A.reshape(dimA, -1).T, B.reshape(dimB, -1).T) + return numpy.allclose(residual, 0) + + +def span_equal(A, B): + # span(A) == span(B) + return span_greater_equal(A, B) and span_greater_equal(B, A) + + +def check_h1div_space(V, degree, reduced=False, bubble=False): + # Test that the divergence of the polynomial space V is spanned by a C0 basis A = V.get_reference_element() top = A.get_topology() sd = A.get_spatial_dimension() + z = (0,)*sd pts = [] for dim in top: for entity in top[dim]: - pts.extend(A.make_points(dim, entity, degree)) + pts.extend(A.make_points(dim, entity, degree+2)) V_tab = V.tabulate(pts, 1) V_div = sum(V_tab[alpha][:, alpha.index(1), :] for alpha in V_tab if sum(alpha) == 1) C0 = CkPolynomialSet(A, degree-1, order=0, variant="bubble") - C0_tab = C0.tabulate(pts)[(0,)*sd] - _, residual, *_ = numpy.linalg.lstsq(C0_tab.T, V_div.T) - assert numpy.allclose(residual, 0) + C0_tab = C0.tabulate(pts)[z] + assert span_equal(V_div, C0_tab) + + if bubble: + # Test that div(Bubbles) = C0 int H^1_0 + assert span_equal(V_div[-(sd+1):], C0_tab[-1:]) + + if not reduced: + # Test that V includes Pk + cell = A.get_parent() or A + Pk = ONPolynomialSet(cell, degree, shape=(sd,)) + Pk_tab = Pk.tabulate(pts)[z] + assert span_greater_equal(V_tab[z], Pk_tab) + + +@pytest.mark.parametrize("cell", (T, S)) +@pytest.mark.parametrize("degree", (2, 3, 4)) +def test_h1div_alfeld_sorokina(cell, degree): + # Test that the divergence of the Alfeld-Sorokina space is spanned by a C0 basis + V = AlfeldSorokinaSpace(cell, degree) + check_h1div_space(V, degree) + + +@pytest.mark.parametrize("cell", (S,)) +@pytest.mark.parametrize("reduced", (False, True), ids=("full", "reduced")) +def test_h1div_guzman_neilan(cell, reduced): + # Test that the divergence of AS + GN Bubble is spanned by a C0 basis + sd = cell.get_spatial_dimension() + degree = 2 + fe = GuzmanNeilanH1div(cell, degree, reduced=reduced) + reduced_dim = fe.space_dimension() - (sd-1)*(sd+1) + V = fe.get_nodal_basis().take(list(range(reduced_dim))) + check_h1div_space(V, degree, reduced=reduced, bubble=True) @pytest.mark.parametrize("cell", (T,)) @@ -77,24 +119,22 @@ def test_hct_stokes_complex(cell, reduced): H1div = sum(H1tab[alpha][:, alpha.index(1), :] for alpha in H1tab if sum(alpha) == 1) H2curl = numpy.stack([H2tab[(0, 1)], -H2tab[(1, 0)]], axis=1) - H2dim = H2.space_dimension() - H1dim = H1.space_dimension() - _, residual, *_ = numpy.linalg.lstsq(H1val.reshape(H1dim, -1).T, H2curl.reshape(H2dim, -1).T) - assert numpy.allclose(residual, 0) - - _, residual, *_ = numpy.linalg.lstsq(L2val.T, H1div.T) - assert numpy.allclose(residual, 0) + assert span_greater_equal(H1val, H2curl) + assert span_greater_equal(L2val, H1div) @pytest.mark.parametrize("cell", (T, S)) -@pytest.mark.parametrize("family", ("AQ", "CH", "GN")) +@pytest.mark.parametrize("family", ("AQ", "CH", "GN", "GN2")) def test_minimal_stokes_space(cell, family): # Test that the C0 Stokes space is spanned by a C0 basis # Also test that its divergence is constant sd = cell.get_spatial_dimension() if family == "GN": degree = 1 - space = ExtendedGuzmanNeilanSpace + space = GuzmanNeilanSpace + elif family == "GN2": + degree = 1 + space = lambda *args, **kwargs: GuzmanNeilanSpace(*args, kind=2, **kwargs) elif family == "CH": degree = 1 space = ChristiansenHuSpace @@ -135,10 +175,14 @@ def test_minimal_stokes_space(cell, family): # Test that divergence is in P0 div = sum(tab[alpha][:, alpha.index(1), :] for alpha in tab if sum(alpha) == 1)[:Vdim] - assert numpy.allclose(div, div[:, 0][:, None]) + if family == "GN2": + # Test that div(GN2) is in P0(Alfeld) + P0 = ONPolynomialSet(K, 0) + P0_tab = P0.tabulate(pts)[(0,)*sd] + assert span_equal(div, P0_tab) + else: + assert numpy.allclose(div, div[:, 0][:, None]) # Test that the full space includes the reduced space assert Wdim > Vdim - _, residual, *_ = numpy.linalg.lstsq(Wtab[z].reshape(Wdim, -1).T, - Vtab[z].reshape(Vdim, -1).T) - assert numpy.allclose(residual, 0) + assert span_greater_equal(Wtab[z], Vtab[z]) From a3a087cb7b4b2e6266a587490d0b5bf9815bb5d9 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 19 Oct 2024 00:14:03 +0100 Subject: [PATCH 07/13] Fix quadratic Guzman-Neilan (#91) * More robust testing for Stokes macroelements --- FIAT/bernardi_raugel.py | 9 +- FIAT/dual_set.py | 8 +- FIAT/guzman_neilan.py | 30 ++--- FIAT/hct.py | 3 +- test/unit/test_stokes_complex.py | 182 ++++++++++++++++++++----------- 5 files changed, 151 insertions(+), 81 deletions(-) diff --git a/FIAT/bernardi_raugel.py b/FIAT/bernardi_raugel.py index 6809dfdd0..c7f4c94ed 100644 --- a/FIAT/bernardi_raugel.py +++ b/FIAT/bernardi_raugel.py @@ -28,12 +28,19 @@ def BernardiRaugelSpace(ref_el, order): entity_ids = expansions.polynomial_entity_ids(ref_el, sd, continuity="C0") slices = {dim: slice(math.comb(order-1, dim)) for dim in range(order)} - slices[sd-1] = slice(None) + slices.pop(sd-1, None) ids = [i + j * dimPd for dim in slices for f in sorted(entity_ids[dim]) for i in entity_ids[dim][f][slices[dim]] for j in range(sd)] + + interior_facets = ref_el.get_interior_facets(sd-1) or () + facets = list(set(entity_ids[sd-1]) - set(interior_facets)) + ids.extend(i + j*dimPd + for f in sorted(facets) + for i in entity_ids[sd-1][f] + for j in range(sd)) return Pd.take(ids) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index 438ffb9af..74630363d 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -115,18 +115,22 @@ def to_riesz(self, poly_set): dpts = set() Qs_to_ells = dict() for i, ell in enumerate(self.nodes): + if len(ell.deriv_dict) > 0: + dpts.update(ell.deriv_dict.keys()) + continue if isinstance(ell, functional.IntegralMoment): Q = ell.Q else: Q = None pts.update(ell.pt_dict.keys()) - dpts.update(ell.deriv_dict.keys()) if Q in Qs_to_ells: Qs_to_ells[Q].append(i) else: Qs_to_ells[Q] = [i] - Qs_to_pts = {None: tuple(sorted(pts))} + Qs_to_pts = {} + if len(pts) > 0: + Qs_to_pts[None] = tuple(sorted(pts)) for Q in Qs_to_ells: if Q is not None: cur_pts = tuple(map(tuple, Q.pts)) diff --git a/FIAT/guzman_neilan.py b/FIAT/guzman_neilan.py index 0c721a308..b11550ed9 100644 --- a/FIAT/guzman_neilan.py +++ b/FIAT/guzman_neilan.py @@ -44,17 +44,17 @@ def GuzmanNeilanSpace(ref_el, order, kind=1, reduced=False): if sd > 2: B = modified_bubble_subspace(B) + K = ref_complex if kind == 2 else ref_el + num_bubbles = sd + 1 if reduced: - BR = BernardiRaugel(ref_el, order).get_nodal_basis() + BR = BernardiRaugel(K, order).get_nodal_basis() reduced_dim = BR.get_num_members() - (sd-1) * (sd+1) BR = BR.take(list(range(reduced_dim))) else: - BR = BernardiRaugelSpace(ref_el, order) + num_bubbles *= sd + BR = BernardiRaugelSpace(K, order) - GN = constant_div_projection(BR, C0, B) - if kind == 2: - Bk = take_interior_bubbles(C0, order) - GN = polynomial_set.polynomial_set_union_normalized(GN, Bk) + GN = constant_div_projection(BR, C0, B, num_bubbles) return GN @@ -108,13 +108,18 @@ def GuzmanNeilanH1div(ref_el, degree=2, reduced=False): This element belongs to a Stokes complex, and is paired with CG1(Alfeld). """ + order = 0 AS = AlfeldSorokina(ref_el, 2) if reduced: - AS = RestrictedElement(AS, restriction_domain="vertex") + order = 1 + # Only extract the div bubbles + div_nodes = [i for i, node in enumerate(AS.dual_basis()) + if len(node.deriv_dict) > 0] + AS = RestrictedElement(AS, indices=div_nodes) elif ref_el.get_spatial_dimension() <= 2: # Quadratic bubbles are already included in 2D return AS - GN = GuzmanNeilanH1(ref_el, order=0) + GN = GuzmanNeilanH1(ref_el, order=order) return NodalEnrichedElement(AS, GN) @@ -189,9 +194,8 @@ def modified_bubble_subspace(B): return M -def constant_div_projection(BR, C0, M): - """Project the BR space into C0 Pk(Alfeld)^d with P_{k-1} divergence. - """ +def constant_div_projection(BR, C0, M, num_bubbles): + """Project the BR space into C0 Pk(Alfeld)^d with P_{k-1} divergence.""" ref_complex = C0.get_reference_element() sd = ref_complex.get_spatial_dimension() degree = C0.degree @@ -208,14 +212,14 @@ def constant_div_projection(BR, C0, M): X = BR.tabulate(qpts, 1) # Invert the divergence B = inner(P, div(U), qwts) - g = inner(P, div(X), qwts) + g = inner(P, div(X)[-num_bubbles:], qwts) w = numpy.linalg.solve(B, g) # Add correction to BR bubbles v = C0.tabulate(qpts)[(0,)*sd] coeffs = numpy.linalg.solve(inner(v, v, qwts), inner(v, X[(0,)*sd], qwts)) coeffs = coeffs.T.reshape(BR.get_num_members(), sd, -1) - coeffs -= numpy.tensordot(w, M.get_coeffs(), axes=(0, 0)) + coeffs[-num_bubbles:] -= numpy.tensordot(w, M.get_coeffs(), axes=(0, 0)) GN = polynomial_set.PolynomialSet(ref_complex, degree, degree, C0.get_expansion_set(), coeffs) return GN diff --git a/FIAT/hct.py b/FIAT/hct.py index b09741cdc..c7c8d6e20 100644 --- a/FIAT/hct.py +++ b/FIAT/hct.py @@ -84,4 +84,5 @@ def __init__(self, ref_el, degree=3, reduced=False): ref_complex = macro.AlfeldSplit(ref_el) dual = HCTDualSet(ref_complex, degree, reduced=reduced) poly_set = macro.CkPolynomialSet(ref_complex, degree, order=1, vorder=degree-1, variant="bubble") - super().__init__(poly_set, dual, degree) + formdegree = 0 + super().__init__(poly_set, dual, degree, formdegree=formdegree) diff --git a/test/unit/test_stokes_complex.py b/test/unit/test_stokes_complex.py index d163ae60c..4abb90322 100644 --- a/test/unit/test_stokes_complex.py +++ b/test/unit/test_stokes_complex.py @@ -13,7 +13,11 @@ from FIAT.alfeld_sorokina import AlfeldSorokinaSpace from FIAT.arnold_qin import ArnoldQinSpace from FIAT.christiansen_hu import ChristiansenHuSpace -from FIAT.guzman_neilan import GuzmanNeilanSpace, GuzmanNeilanH1div +from FIAT.guzman_neilan import ( + GuzmanNeilanFirstKindH1, + GuzmanNeilanSecondKindH1, + GuzmanNeilanH1div, + GuzmanNeilanSpace) from FIAT.restricted import RestrictedElement @@ -21,11 +25,20 @@ S = ufc_simplex(3) +def rHCT(cell): + return RestrictedElement(HCT(cell, reduced=True), restriction_domain="vertex") + + +def rAQ(cell): + return RestrictedElement(AQ(cell, reduced=True), indices=list(range(9))) + + def span_greater_equal(A, B): # span(A) >= span(B) dimA = A.shape[0] dimB = B.shape[0] - _, residual, *_ = numpy.linalg.lstsq(A.reshape(dimA, -1).T, B.reshape(dimB, -1).T) + _, residual, *_ = numpy.linalg.lstsq(A.reshape(dimA, -1).T, + B.reshape(dimB, -1).T) return numpy.allclose(residual, 0) @@ -34,21 +47,34 @@ def span_equal(A, B): return span_greater_equal(A, B) and span_greater_equal(B, A) +def div(U): + """Return divergence from dict of tabulations """ + return sum(U[k][:, k.index(1), :] for k in U if sum(k) == 1) + + +def rot(U): + """Return rot from dict of tabulations """ + return numpy.stack([U[(0, 1)], -U[(1, 0)]], axis=1) + + +def make_points(K, degree): + top = K.get_topology() + pts = [] + for dim in top: + for entity in top[dim]: + pts.extend(K.make_points(dim, entity, degree)) + return pts + + def check_h1div_space(V, degree, reduced=False, bubble=False): # Test that the divergence of the polynomial space V is spanned by a C0 basis A = V.get_reference_element() - top = A.get_topology() sd = A.get_spatial_dimension() z = (0,)*sd - pts = [] - for dim in top: - for entity in top[dim]: - pts.extend(A.make_points(dim, entity, degree+2)) - + pts = make_points(A, degree+2) V_tab = V.tabulate(pts, 1) - V_div = sum(V_tab[alpha][:, alpha.index(1), :] - for alpha in V_tab if sum(alpha) == 1) + V_div = div(V_tab) C0 = CkPolynomialSet(A, degree-1, order=0, variant="bubble") C0_tab = C0.tabulate(pts)[z] @@ -58,12 +84,12 @@ def check_h1div_space(V, degree, reduced=False, bubble=False): # Test that div(Bubbles) = C0 int H^1_0 assert span_equal(V_div[-(sd+1):], C0_tab[-1:]) - if not reduced: - # Test that V includes Pk - cell = A.get_parent() or A - Pk = ONPolynomialSet(cell, degree, shape=(sd,)) - Pk_tab = Pk.tabulate(pts)[z] - assert span_greater_equal(V_tab[z], Pk_tab) + k = degree - 1 if reduced else degree + # Test that V includes Pk + cell = A.get_parent() or A + Pk = ONPolynomialSet(cell, k, shape=(sd,)) + Pk_tab = Pk.tabulate(pts)[z] + assert span_greater_equal(V_tab[z], Pk_tab) @pytest.mark.parametrize("cell", (T, S)) @@ -86,41 +112,74 @@ def test_h1div_guzman_neilan(cell, reduced): check_h1div_space(V, degree, reduced=reduced, bubble=True) -@pytest.mark.parametrize("cell", (T,)) -@pytest.mark.parametrize("reduced", (False, True)) -def test_hct_stokes_complex(cell, reduced): - # Test that we have the lowest-order discrete Stokes complex, verifying - # that the range of the exterior derivative of each space is contained by - # the next space in the sequence - H2 = HCT(T, reduced=reduced) - A = H2.get_reference_complex() - if reduced: - # HCT-red(3) -curl-> AQ-red(2) -div-> DG(0) - H2 = RestrictedElement(H2, restriction_domain="vertex") - H1 = RestrictedElement(AQ(T, reduced=reduced), indices=list(range(9))) - L2 = DG(cell, 0) - else: - # HCT(3) -curl-> AS(2) -div-> CG(1, Alfeld) - H1 = AS(cell) - L2 = CG(A, 1) +def check_stokes_complex(spaces, degree): + # Test that we have a discrete Stokes complex, verifying that the range of + # the exterior derivative of each space is contained by the next space in + # the sequence + A = spaces[0].get_reference_complex() + sd = A.get_spatial_dimension() + z = (0,) * sd + + pts = make_points(A, degree+2) + tab = [V.tabulate(1, pts) for V in spaces] + if len(tab) > 2: + # check rot(V0) in V1 + assert span_greater_equal(tab[1][z], rot(tab[0])) + + # check div(V1) = V2 + assert span_equal(tab[-1][z], div(tab[-2])) + + # Test that V1 includes Pk + cell = A.get_parent() or A + Pk = ONPolynomialSet(cell, degree, shape=(sd,)) + Pk_tab = Pk.tabulate(pts)[z] + assert span_greater_equal(tab[-2][z], Pk_tab) - pts = [] - top = A.get_topology() - for dim in top: - for entity in top[dim]: - pts.extend(A.make_points(dim, entity, 4)) - H2tab = H2.tabulate(1, pts) - H1tab = H1.tabulate(1, pts) - L2tab = L2.tabulate(0, pts) +@pytest.mark.parametrize("reduced", (False, True), ids=("full", "reduced")) +@pytest.mark.parametrize("sobolev", ("H1", "H1div")) +@pytest.mark.parametrize("cell", (T,)) +def test_hct_stokes_complex(cell, sobolev, reduced): + if sobolev == "H1": + if reduced: + spaces = [rHCT(cell), rAQ(cell), DG(cell, 0)] + degree = 1 + else: + spaces = [HCT(cell), AQ(cell), DG(cell, 0)] + degree = 1 + elif sobolev == "H1div": + if reduced: + spaces = [rHCT(cell), + GuzmanNeilanH1div(cell, reduced=True), + CG(cell, 1, variant="alfeld")] + degree = 1 + else: + spaces = [HCT(cell), AS(cell), CG(cell, 1, variant="alfeld")] + degree = 2 + else: + raise ValueError(f"Unexpected sobolev space {sobolev}") + check_stokes_complex(spaces, degree) - L2val = L2tab[(0, 0)] - H1val = H1tab[(0, 0)] - H1div = sum(H1tab[alpha][:, alpha.index(1), :] for alpha in H1tab if sum(alpha) == 1) - H2curl = numpy.stack([H2tab[(0, 1)], -H2tab[(1, 0)]], axis=1) - assert span_greater_equal(H1val, H2curl) - assert span_greater_equal(L2val, H1div) +@pytest.mark.parametrize("cell", (T, S)) +@pytest.mark.parametrize("kind", (1, 2, "H1div", "H1div-red")) +def test_gn_stokes_pairs(cell, kind): + order = cell.get_spatial_dimension() - 1 + if kind == 1: + spaces = [GuzmanNeilanFirstKindH1(cell, order), DG(cell, order-1)] + degree = order + elif kind == 2: + spaces = [GuzmanNeilanSecondKindH1(cell, order), DG(cell, order-1, variant="alfeld")] + degree = order + elif kind == "H1div": + spaces = [GuzmanNeilanH1div(cell), CG(cell, 1, variant="alfeld")] + degree = 2 + elif kind == "H1div-red": + spaces = [GuzmanNeilanH1div(cell, reduced=True), CG(cell, 1, variant="alfeld")] + degree = 1 + else: + raise ValueError(f"Unexpected kind {kind}") + check_stokes_complex(spaces, degree) @pytest.mark.parametrize("cell", (T, S)) @@ -150,38 +209,33 @@ def test_minimal_stokes_space(cell, family): Vdim = V.get_num_members() K = W.get_reference_element() sd = K.get_spatial_dimension() - top = K.get_topology() - - pts = [] - for dim in top: - for entity in top[dim]: - pts.extend(K.make_points(dim, entity, degree)) + pts = make_points(K, degree+2) C0 = CkPolynomialSet(K, sd, order=0, variant="bubble") C0_tab = C0.tabulate(pts) Wtab = W.tabulate(pts, 1) Vtab = V.tabulate(pts, 1) - z = (0,)*sd - for tab in (Vtab, Wtab): + z = (0,) * sd + for Xtab in (Vtab, Wtab): # Test that the space is full rank - _, sig, _ = numpy.linalg.svd(tab[z].reshape(-1, sd*len(pts)).T, full_matrices=True) + _, sig, _ = numpy.linalg.svd(Xtab[z].reshape(-1, sd*len(pts)).T, full_matrices=True) assert all(sig > 1E-10) # Test that the space is C0 for k in range(sd): - _, residual, *_ = numpy.linalg.lstsq(C0_tab[z].T, tab[z][:, k, :].T) + _, residual, *_ = numpy.linalg.lstsq(C0_tab[z].T, Xtab[z][:, k, :].T) assert numpy.allclose(residual, 0) # Test that divergence is in P0 - div = sum(tab[alpha][:, alpha.index(1), :] - for alpha in tab if sum(alpha) == 1)[:Vdim] - if family == "GN2": + divX = div(Xtab)[:Vdim] + if family in {"GN", "GN2"}: # Test that div(GN2) is in P0(Alfeld) - P0 = ONPolynomialSet(K, 0) - P0_tab = P0.tabulate(pts)[(0,)*sd] - assert span_equal(div, P0_tab) + ref_el = K if family == "GN2" else K.get_parent() + P0 = ONPolynomialSet(ref_el, degree-1) + P0_tab = P0.tabulate(pts)[z] + assert span_equal(divX, P0_tab) else: - assert numpy.allclose(div, div[:, 0][:, None]) + assert numpy.allclose(divX, divX[:, 0][:, None]) # Test that the full space includes the reduced space assert Wdim > Vdim From 1088bfd666261c7043e5fde3bbebb5ffd87ac0bd Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 20 Oct 2024 14:21:58 +0100 Subject: [PATCH 08/13] Cleanup Regge and HHJ (#92) --- FIAT/functional.py | 37 ++++++----- FIAT/hellan_herrmann_johnson.py | 107 +++++++++++++----------------- FIAT/regge.py | 112 ++++++++++---------------------- 3 files changed, 101 insertions(+), 155 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index 36989e607..1284cd5e2 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -170,7 +170,7 @@ class PointEvaluation(Functional): def __init__(self, ref_el, x): pt_dict = {x: [(1.0, tuple())]} - Functional.__init__(self, ref_el, tuple(), pt_dict, {}, "PointEval") + super().__init__(ref_el, tuple(), pt_dict, {}, "PointEval") def __call__(self, fn): """Evaluate the functional on the function fn.""" @@ -183,17 +183,18 @@ def tostr(self): class ComponentPointEvaluation(Functional): """Class representing point evaluation of a particular component - of a vector function at a particular point x.""" + of a vector/tensor function at a particular point x.""" def __init__(self, ref_el, comp, shp, x): - if len(shp) != 1: - raise Exception("Illegal shape") - if comp < 0 or comp >= shp[0]: - raise Exception("Illegal component") + if not isinstance(comp, tuple): + comp = (comp,) + if len(shp) != len(comp): + raise ValueError("Component and shape are incompatible") + if any(i < 0 or i >= n for i, n in zip(comp, shp)): + raise ValueError("Illegal component") self.comp = comp - pt_dict = {x: [(1.0, (comp,))]} - Functional.__init__(self, ref_el, shp, pt_dict, {}, - "ComponentPointEval") + pt_dict = {x: [(1.0, comp)]} + super().__init__(ref_el, shp, pt_dict, {}, "ComponentPointEval") def tostr(self): x = list(map(str, list(self.pt_dict.keys())[0])) @@ -209,7 +210,7 @@ def __init__(self, ref_el, x, alpha): self.alpha = tuple(alpha) self.order = sum(self.alpha) - Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointDeriv") + super().__init__(ref_el, tuple(), {}, dpt_dict, "PointDeriv") def __call__(self, fn): """Evaluate the functional on the function fn. Note that this depends @@ -239,7 +240,7 @@ def __init__(self, ref_el, facet_no, pt): alphas.append(alpha) dpt_dict = {pt: [(n[i], tuple(alphas[i]), tuple()) for i in range(sd)]} - Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv") + super().__init__(ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv") class PointNormalSecondDerivative(Functional): @@ -265,7 +266,7 @@ def __init__(self, ref_el, facet_no, pt): self.alphas = alphas dpt_dict = {pt: [(n[i], alphas[i], tuple()) for i in range(sd)]} - Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv") + super().__init__(ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv") class PointDivergence(Functional): @@ -273,9 +274,11 @@ class PointDivergence(Functional): functions at a particular point x.""" def __init__(self, ref_el, x): - dpt_dict = {x: [(1.0, alpha, (alpha.index(1),)) for alpha in polynomial_set.mis(len(x), 1)]} + sd = ref_el.get_spatial_dimension() + alphas = tuple(map(tuple, numpy.eye(sd, dtype=int))) + dpt_dict = {x: [(1.0, alpha, (alpha.index(1),)) for alpha in alphas]} - Functional.__init__(self, ref_el, (len(x),), {}, dpt_dict, "PointDiv") + super().__init__(ref_el, (len(x),), {}, dpt_dict, "PointDiv") class IntegralMoment(Functional): @@ -297,7 +300,7 @@ def __init__(self, ref_el, Q, f_at_qpts, comp=tuple(), shp=tuple()): self.comp = comp weights = numpy.multiply(f_at_qpts, qwts) pt_dict = {tuple(pt): [(wt, comp)] for pt, wt in zip(qpts, weights)} - Functional.__init__(self, ref_el, shp, pt_dict, {}, "IntegralMoment") + super().__init__(ref_el, shp, pt_dict, {}, "IntegralMoment") def __call__(self, fn): """Evaluate the functional on the function fn.""" @@ -331,8 +334,8 @@ def __init__(self, ref_el, facet_no, Q, f_at_qpts): dpt_dict = {tuple(pt): [(wt*n[i], alphas[i], tuple()) for i in range(sd)] for pt, wt in zip(points, weights)} - Functional.__init__(self, ref_el, tuple(), - {}, dpt_dict, "IntegralMomentOfNormalDerivative") + super().__init__(ref_el, tuple(), + {}, dpt_dict, "IntegralMomentOfNormalDerivative") class IntegralLegendreDirectionalMoment(Functional): diff --git a/FIAT/hellan_herrmann_johnson.py b/FIAT/hellan_herrmann_johnson.py index 8218b00f6..631272c64 100644 --- a/FIAT/hellan_herrmann_johnson.py +++ b/FIAT/hellan_herrmann_johnson.py @@ -7,90 +7,75 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT.finite_element import CiarletElement -from FIAT.dual_set import DualSet -from FIAT.polynomial_set import ONSymTensorPolynomialSet -from FIAT.functional import PointwiseInnerProductEvaluation as InnerProduct -import numpy +from FIAT import dual_set, finite_element, polynomial_set +from FIAT.functional import ComponentPointEvaluation, PointwiseInnerProductEvaluation as InnerProduct -class HellanHerrmannJohnsonDual(DualSet): +class HellanHerrmannJohnsonDual(dual_set.DualSet): """Degrees of freedom for Hellan-Herrmann-Johnson elements.""" - def __init__(self, cell, degree): - dim = cell.get_spatial_dimension() - if not dim == 2: - raise ValueError("Hellan_Herrmann-Johnson elements are only" + def __init__(self, ref_el, degree): + sd = ref_el.get_spatial_dimension() + if sd != 2: + raise ValueError("Hellan-Herrmann-Johnson elements are only" "defined in dimension 2.") - # construct the degrees of freedoms - dofs = [] # list of functionals - # dof_ids[i][j] contains the indices of dofs that are associated with - # entity j in dim i - dof_ids = {} + top = ref_el.get_topology() + entity_ids = {dim: {i: [] for i in sorted(top[dim])} for dim in sorted(top)} + nodes = [] - # no vertex dof - dof_ids[0] = {i: [] for i in range(dim + 1)} - # edge dofs - (_dofs, _dof_ids) = self._generate_edge_dofs(cell, degree, 0) - dofs.extend(_dofs) - dof_ids[1] = _dof_ids - # cell dofs - (_dofs, _dof_ids) = self._generate_trig_dofs(cell, degree, len(dofs)) - dofs.extend(_dofs) - dof_ids[dim] = _dof_ids + # no vertex dofs + edge_dofs, entity_ids[1] = self._generate_edge_dofs(ref_el, degree, 0) + nodes.extend(edge_dofs) + cell_nodes, entity_ids[sd] = self._generate_cell_dofs(ref_el, degree, len(nodes)) + nodes.extend(cell_nodes) - super().__init__(dofs, cell, dof_ids) + super().__init__(nodes, ref_el, entity_ids) @staticmethod - def _generate_edge_dofs(cell, degree, offset): + def _generate_edge_dofs(ref_el, degree, offset): """Generate dofs on edges. On each edge, let n be its normal. For degree=r, the scalar function n^T u n - is evaluated at points enough to control P(r). + is evaluated at enough points to control P(r). """ - dofs = [] - dof_ids = {} - for entity_id in range(3): # a triangle has 3 edges - pts = cell.make_points(1, entity_id, degree + 2) # edges are 1D - normal = cell.compute_scaled_normal(entity_id) - dofs += [InnerProduct(cell, normal, normal, pt) for pt in pts] - num_new_dofs = len(pts) # 1 dof per point on edge - dof_ids[entity_id] = list(range(offset, offset + num_new_dofs)) - offset += num_new_dofs - return (dofs, dof_ids) + dim = 1 + top = ref_el.get_topology() + nodes = [] + entity_ids = {} + for entity_id in sorted(top[dim]): + pts = ref_el.make_points(dim, entity_id, degree + 2) + normal = ref_el.compute_scaled_normal(entity_id) + nodes.extend(InnerProduct(ref_el, normal, normal, pt) for pt in pts) + num_new_nodes = len(pts) + entity_ids[entity_id] = list(range(offset, offset + num_new_nodes)) + offset += num_new_nodes + return nodes, entity_ids @staticmethod - def _generate_trig_dofs(cell, degree, offset): - """Generate dofs on edges. + def _generate_cell_dofs(ref_el, degree, offset): + """Generate dofs on the cell interior. On each triangle, for degree=r, the three components u11, u12, u22 - are evaluated at points enough to control P(r-1). + are evaluated at enough points to control P(r-1). """ - dofs = [] - dof_ids = {} - pts = cell.make_points(2, 0, degree + 2) # 2D trig #0 - e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 - e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 - basis = [(e1, e1), (e1, e2), (e2, e2)] # basis for symmetric matrix - for (v1, v2) in basis: - dofs += [InnerProduct(cell, v1, v2, pt) for pt in pts] - num_dofs = 3 * len(pts) # 3 dofs per trig - dof_ids[0] = list(range(offset, offset + num_dofs)) - return (dofs, dof_ids) + sd = ref_el.get_spatial_dimension() + shp = (sd, sd) + basis = [(i, j) for i in range(sd) for j in range(i, sd)] + pts = ref_el.make_points(sd, 0, degree + 2) + nodes = [ComponentPointEvaluation(ref_el, comp, shp, pt) + for comp in basis for pt in pts] + entity_ids = {0: list(range(offset, offset + len(nodes)))} + return nodes, entity_ids -class HellanHerrmannJohnson(CiarletElement): +class HellanHerrmannJohnson(finite_element.CiarletElement): """The definition of Hellan-Herrmann-Johnson element. It is defined only in dimension 2. It consists of piecewise polynomial symmetric-matrix-valued functions of degree r or less with normal-normal continuity. """ - def __init__(self, cell, degree): + def __init__(self, ref_el, degree): assert degree >= 0, "Hellan-Herrmann-Johnson starts at degree 0!" - # shape functions - Ps = ONSymTensorPolynomialSet(cell, degree) - # degrees of freedom - Ls = HellanHerrmannJohnsonDual(cell, degree) - # mapping under affine transformation + poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + dual = HellanHerrmannJohnsonDual(ref_el, degree) mapping = "double contravariant piola" - - super().__init__(Ps, Ls, degree, mapping=mapping) + super().__init__(poly_set, dual, degree, mapping=mapping) diff --git a/FIAT/regge.py b/FIAT/regge.py index e9292f663..6e08d0399 100644 --- a/FIAT/regge.py +++ b/FIAT/regge.py @@ -6,90 +6,48 @@ # This file is part of FIAT (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT.finite_element import CiarletElement -from FIAT.dual_set import DualSet -from FIAT.polynomial_set import ONSymTensorPolynomialSet +from FIAT import dual_set, finite_element, polynomial_set from FIAT.functional import PointwiseInnerProductEvaluation as InnerProduct -class ReggeDual(DualSet): - """Degrees of freedom for generalized Regge finite elements.""" - def __init__(self, cell, degree): - dim = cell.get_spatial_dimension() - if (dim < 2) or (dim > 3): - raise ValueError("Generalized Regge elements are implemented only " - "for dimension 2--3. For 1D, it is just DG(r).") +class ReggeDual(dual_set.DualSet): + """Degrees of freedom for generalized Regge finite elements. - # construct the degrees of freedoms - dofs = [] # list of functionals - # dof_ids[i][j] contains the indices of dofs that are associated with - # entity j in dim i - dof_ids = {} - - # no vertex dof - dof_ids[0] = {i: [] for i in range(dim + 1)} - # edge dofs - (_dofs, _dof_ids) = self._generate_dofs(cell, 1, degree, 0) - dofs.extend(_dofs) - dof_ids[1] = _dof_ids - # facet dofs for 3D - if dim == 3: - (_dofs, _dof_ids) = self._generate_dofs(cell, 2, degree, len(dofs)) - dofs.extend(_dofs) - dof_ids[2] = _dof_ids - # cell dofs - (_dofs, _dof_ids) = self._generate_dofs(cell, dim, degree, len(dofs)) - dofs.extend(_dofs) - dof_ids[dim] = _dof_ids - - super().__init__(dofs, cell, dof_ids) - - @staticmethod - def _generate_dofs(cell, entity_dim, degree, offset): - """Generate degrees of freedom for enetities of dimension entity_dim - - Input: all obvious except - offset -- the current first available dof id. - - Output: - dofs -- an array of dofs associated to entities in that dim - dof_ids -- a dict mapping entity_id to the range of indices of dofs - associated to it. - - On a k-face for degree r, the dofs are given by the value of - t^T u t - evaluated at points enough to control P(r-k+1) for all the edge - tangents of the face. - `cell.make_points(entity_dim, entity_id, degree + 2)` happens to - generate exactly those points needed. - """ - dofs = [] - dof_ids = {} - num_entities = len(cell.get_topology()[entity_dim]) - for entity_id in range(num_entities): - pts = cell.make_points(entity_dim, entity_id, degree + 2) - tangents = cell.compute_face_edge_tangents(entity_dim, entity_id) - dofs += [InnerProduct(cell, t, t, pt) - for pt in pts - for t in tangents] - num_new_dofs = len(pts) * len(tangents) - dof_ids[entity_id] = list(range(offset, offset + num_new_dofs)) - offset += num_new_dofs - return (dofs, dof_ids) - - -class Regge(CiarletElement): + On a k-face for degree r, the dofs are given by the value of + t^T u t + evaluated at enough points to control P(r-k+1) for all the edge + tangents of the face. + `ref_el.make_points(dim, entity, degree + 2)` happens to + generate exactly those points needed. + """ + def __init__(self, ref_el, degree): + top = ref_el.get_topology() + entity_ids = {dim: {i: [] for i in sorted(top[dim])} for dim in sorted(top)} + nodes = [] + for dim in sorted(top): + if dim == 0: + # no vertex dofs + continue + for entity in sorted(top[dim]): + cur = len(nodes) + tangents = ref_el.compute_face_edge_tangents(dim, entity) + pts = ref_el.make_points(dim, entity, degree + 2) + nodes.extend(InnerProduct(ref_el, t, t, pt) + for pt in pts + for t in tangents) + entity_ids[dim][entity].extend(range(cur, len(nodes))) + + super().__init__(nodes, ref_el, entity_ids) + + +class Regge(finite_element.CiarletElement): """The generalized Regge elements for symmetric-matrix-valued functions. REG(r) in dimension n is the space of polynomial symmetric-matrix-valued functions of degree r or less with tangential-tangential continuity. """ - def __init__(self, cell, degree): + def __init__(self, ref_el, degree): assert degree >= 0, "Regge start at degree 0!" - # shape functions - Ps = ONSymTensorPolynomialSet(cell, degree) - # degrees of freedom - Ls = ReggeDual(cell, degree) - # mapping under affine transformation + poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + dual = ReggeDual(ref_el, degree) mapping = "double covariant piola" - - super().__init__(Ps, Ls, degree, mapping=mapping) + super().__init__(poly_set, dual, degree, mapping=mapping) From 49e6bdc1bdb3dfa41c8c91c72b2e8832e01f4a3c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 23 Oct 2024 09:32:15 +0100 Subject: [PATCH 09/13] more robust tests for HCT/HCT-red (#94) --- test/unit/test_hct.py | 47 ++++++++++++++++++++++++++++++-- test/unit/test_stokes_complex.py | 6 ++-- 2 files changed, 46 insertions(+), 7 deletions(-) diff --git a/test/unit/test_hct.py b/test/unit/test_hct.py index 4dbb96024..84e0a2379 100644 --- a/test/unit/test_hct.py +++ b/test/unit/test_hct.py @@ -1,9 +1,10 @@ import pytest import numpy -from FIAT import HsiehCloughTocher as HCT -from FIAT.reference_element import ufc_simplex, make_lattice +from FIAT import RestrictedElement, HsiehCloughTocher as HCT +from FIAT.reference_element import ufc_simplex from FIAT.functional import PointEvaluation +from FIAT.macro import CkPolynomialSet @pytest.fixture @@ -13,12 +14,28 @@ def cell(): return K +def span_greater_equal(A, B): + # span(A) >= span(B) + _, residual, *_ = numpy.linalg.lstsq(A.reshape(A.shape[0], -1).T, + B.reshape(B.shape[0], -1).T) + return numpy.allclose(residual, 0) + + +def make_points(K, degree): + top = K.get_topology() + pts = [] + for dim in top: + for entity in top[dim]: + pts.extend(K.make_points(dim, entity, degree)) + return pts + + @pytest.mark.parametrize("reduced", (False, True)) def test_hct_constant(cell, reduced): # Test that bfs associated with point evaluation sum up to 1 fe = HCT(cell, reduced=reduced) - pts = make_lattice(cell.get_vertices(), 3) + pts = make_points(cell, 4) tab = fe.tabulate(2, pts) coefs = numpy.zeros((fe.space_dimension(),)) @@ -33,3 +50,27 @@ def test_hct_constant(cell, reduced): expected = 1 if sum(alpha) == 0 else 0 vals = numpy.dot(coefs, tab[alpha]) assert numpy.allclose(vals, expected) + + +@pytest.mark.parametrize("reduced", (False, True)) +def test_full_polynomials(cell, reduced): + # Test that HCT/HCT-red contains all cubics/quadratics + fe = HCT(cell, reduced=reduced) + if reduced: + fe = RestrictedElement(fe, restriction_domain="vertex") + + ref_complex = fe.get_reference_complex() + pts = make_points(ref_complex, 4) + tab = fe.tabulate(0, pts)[(0, 0)] + + degree = fe.degree() + if reduced: + degree -= 1 + + P = CkPolynomialSet(cell, degree, variant="bubble") + P_tab = P.tabulate(pts)[(0, 0)] + assert span_greater_equal(tab, P_tab) + + C1 = CkPolynomialSet(ref_complex, degree, order=1, variant="bubble") + C1_tab = C1.tabulate(pts)[(0, 0)] + assert span_greater_equal(tab, C1_tab) diff --git a/test/unit/test_stokes_complex.py b/test/unit/test_stokes_complex.py index 4abb90322..a3bf36370 100644 --- a/test/unit/test_stokes_complex.py +++ b/test/unit/test_stokes_complex.py @@ -35,10 +35,8 @@ def rAQ(cell): def span_greater_equal(A, B): # span(A) >= span(B) - dimA = A.shape[0] - dimB = B.shape[0] - _, residual, *_ = numpy.linalg.lstsq(A.reshape(dimA, -1).T, - B.reshape(dimB, -1).T) + _, residual, *_ = numpy.linalg.lstsq(A.reshape(A.shape[0], -1).T, + B.reshape(B.shape[0], -1).T) return numpy.allclose(residual, 0) From ae368208aa000843e656b957295d0739727e448c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 25 Oct 2024 20:44:32 +0100 Subject: [PATCH 10/13] Fix GN modified bubble space (#97) --- FIAT/guzman_neilan.py | 21 ++++++++------------- 1 file changed, 8 insertions(+), 13 deletions(-) diff --git a/FIAT/guzman_neilan.py b/FIAT/guzman_neilan.py index b11550ed9..ff853e583 100644 --- a/FIAT/guzman_neilan.py +++ b/FIAT/guzman_neilan.py @@ -17,7 +17,6 @@ from FIAT.quadrature_schemes import create_quadrature from FIAT.restricted import RestrictedElement from FIAT.nodal_enriched import NodalEnrichedElement -from itertools import chain import numpy import math @@ -168,21 +167,17 @@ def modified_bubble_subspace(B): hat = B.take([0]) hat_at_qpts = hat.tabulate(qpts)[(0,)*sd][0, 0] - # tabulate the BDM facet functions - ref_el = ref_complex.get_parent() - BDM = BrezziDouglasMarini(ref_el, degree-1) - entity_dofs = BDM.entity_dofs() - facet_dofs = list(range(BDM.space_dimension() - len(entity_dofs[sd][0]))) - BDM_facet = BDM.get_nodal_basis().take(facet_dofs) - phis = BDM_facet.tabulate(qpts)[(0,)*sd] - # tabulate the bubbles = hat ** (degree - k) * BDMk_facet + ref_el = ref_complex.get_parent() bubbles = [numpy.eye(sd)[:, :, None] * hat_at_qpts[None, None, :] ** degree] for k in range(1, degree): - dimPk = math.comb(k + sd-1, sd-1) - idsPk = list(chain.from_iterable(entity_dofs[sd-1][f][:dimPk] - for f in entity_dofs[sd-1])) - bubbles.append(numpy.multiply(phis[idsPk], hat_at_qpts ** (degree-k))) + # tabulate the BDM facet functions + BDM = BrezziDouglasMarini(ref_el, k) + BDM_facet = BDM.get_nodal_basis().take(BDM.dual.get_indices("facet")) + phis = BDM_facet.tabulate(qpts)[(0,)*sd] + + bubbles.append(numpy.multiply(phis, hat_at_qpts ** (degree-k))) + bubbles = numpy.concatenate(bubbles, axis=0) # store the bubbles into a PolynomialSet via L2 projection From bc9c29d4a8ee9aa2cadd5ffceef868ed93f603ae Mon Sep 17 00:00:00 2001 From: ksagiyam Date: Wed, 4 Sep 2024 00:00:20 +0100 Subject: [PATCH 11/13] quadrature: add quad point symmetry info orientation: add Orientation class --- FIAT/orientation_utils.py | 21 +++++ FIAT/quadrature.py | 35 ++++++++ FIAT/reference_element.py | 172 +++++++++++++++++++++++++++++++++++++- 3 files changed, 227 insertions(+), 1 deletion(-) diff --git a/FIAT/orientation_utils.py b/FIAT/orientation_utils.py index 92c678991..2852f33f4 100644 --- a/FIAT/orientation_utils.py +++ b/FIAT/orientation_utils.py @@ -4,6 +4,27 @@ import numpy as np +class Orientation(object): + """Base class representing unsigned integer orientations. + + Orientations represented by this class are to be consistent with those used in + `make_entity_permutations_simplex` and `make_entity_permutations_tensorproduct`. + + """ + + def __floordiv__(self, other): + raise NotImplementedError("Subclass must implement __floordiv__") + + def __rfloordiv__(self, other): + raise NotImplementedError("Subclass must implement __rfloordiv__") + + def __mod__(self, other): + raise NotImplementedError("Subclass must implement __mod__") + + def __rmod__(self, other): + raise NotImplementedError("Subclass must implement __rmod__") + + def make_entity_permutations_simplex(dim, npoints): r"""Make orientation-permutation map for the given simplex dimension, dim, and the number of points along diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py index 61702c54c..51e4f45e2 100644 --- a/FIAT/quadrature.py +++ b/FIAT/quadrature.py @@ -8,10 +8,12 @@ # Modified by David A. Ham (david.ham@imperial.ac.uk), 2015 import itertools +from math import factorial import numpy from recursivenodes.quadrature import gaussjacobi, lobattogaussjacobi, simplexgausslegendre from FIAT import reference_element +from FIAT.orientation_utils import make_entity_permutations_simplex def pseudo_determinant(A): @@ -51,6 +53,7 @@ def __init__(self, ref_el, pts, wts): self.ref_el = ref_el self.pts = pts self.wts = wts + self._intrinsic_orientation_permutation_map_tuple = (None, ) def get_points(self): return numpy.array(self.pts) @@ -61,6 +64,32 @@ def get_weights(self): def integrate(self, f): return sum(w * f(x) for x, w in zip(self.pts, self.wts)) + @property + def extrinsic_orientation_permutation_map(self): + """A map from extrinsic orientations to corresponding axis permutation matrices. + + Notes + ----- + result[eo] gives the physical axis-reference axis permutation matrix corresponding to + eo (extrinsic orientation). + + """ + return self.ref_el.extrinsic_orientation_permutation_map + + @property + def intrinsic_orientation_permutation_map_tuple(self): + """A tuple of maps from intrinsic orientations to corresponding point permutations for each reference cell axis. + + Notes + ----- + result[axis][io] gives the physical point-reference point permutation array corresponding to + io (intrinsic orientation) on ``axis``. + + """ + if any(m is None for m in self._intrinsic_orientation_permutation_map_tuple): + raise ValueError("Must set _intrinsic_orientation_permutation_map_tuple") + return self._intrinsic_orientation_permutation_map_tuple + class GaussJacobiQuadratureLineRule(QuadratureRule): """Gauss-Jacobi quadature rule determined by Jacobi weights a and b @@ -71,6 +100,12 @@ def __init__(self, ref_el, m, a=0, b=0): pts_ref, wts_ref = gaussjacobi(m, a, b) pts, wts = map_quadrature(pts_ref, wts_ref, Ref1, ref_el) QuadratureRule.__init__(self, ref_el, pts, wts) + # Set _intrinsic_orientation_permutation_map_tuple. + dim = 1 + a = numpy.zeros((factorial(dim + 1), m), dtype=int) + for io, perm in make_entity_permutations_simplex(dim, m).items(): + a[io, perm] = range(m) + self._intrinsic_orientation_permutation_map_tuple = (a, ) class GaussLobattoLegendreQuadratureLineRule(QuadratureRule): diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index f09751578..95f0d8f8a 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -29,8 +29,11 @@ from recursivenodes.nodes import _decode_family, _recursive from FIAT.orientation_utils import ( + Orientation, make_cell_orientation_reflection_map_simplex, - make_cell_orientation_reflection_map_tensorproduct) + make_cell_orientation_reflection_map_tensorproduct, + make_entity_permutations_simplex, +) POINT = 0 LINE = 1 @@ -256,6 +259,52 @@ def cell_orientation_reflection_map(self): """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" raise NotImplementedError("Should be implemented in a subclass.") + def extract_extrinsic_orientation(self, o): + """Extract extrinsic orientation. + + Parameters + ---------- + o : Orientation + Total orientation. + + Returns + ------- + Orientation + Extrinsic orientation. + + """ + raise NotImplementedError("Should be implemented in a subclass.") + + def extract_intrinsic_orientation(self, o, axis): + """Extract intrinsic orientation. + + Parameters + ---------- + o : Orientation + Total orientation. + axis : int + Reference cell axis for which intrinsic orientation is computed. + + Returns + ------- + Orientation + Intrinsic orientation. + + """ + raise NotImplementedError("Should be implemented in a subclass.") + + @property + def extrinsic_orientation_permutation_map(self): + """A map from extrinsic orientations to corresponding axis permutation matrices. + + Notes + ----- + result[eo] gives the physical axis-reference axis permutation matrix corresponding to + eo (extrinsic orientation). + + """ + raise NotImplementedError("Should be implemented in a subclass.") + def is_simplex(self): return False @@ -725,6 +774,58 @@ def contains_point(self, point, epsilon=0.0, entity=None): """ return self.distance_to_point_l1(point, entity=entity) <= epsilon + def extract_extrinsic_orientation(self, o): + """Extract extrinsic orientation. + + Parameters + ---------- + o : Orientation + Total orientation. + + Returns + ------- + Orientation + Extrinsic orientation. + + """ + if not isinstance(o, Orientation): + raise TypeError(f"Expecting an instance of Orientation : got {o}") + return 0 + + def extract_intrinsic_orientation(self, o, axis): + """Extract intrinsic orientation. + + Parameters + ---------- + o : Orientation + Total orientation. + axis : int + Reference cell axis for which intrinsic orientation is computed. + + Returns + ------- + Orientation + Intrinsic orientation. + + """ + if not isinstance(o, Orientation): + raise TypeError(f"Expecting an instance of Orientation : got {o}") + if axis != 0: + raise ValueError(f"axis ({axis}) != 0") + return o + + @property + def extrinsic_orientation_permutation_map(self): + """A map from extrinsic orientations to corresponding axis permutation matrices. + + Notes + ----- + result[eo] gives the physical axis-reference axis permutation matrix corresponding to + eo (extrinsic orientation). + + """ + return numpy.diag((1, )).astype(int).reshape((1, 1, 1)) + class Simplex(SimplicialComplex): r"""Abstract class for a reference simplex. @@ -1162,6 +1263,75 @@ def __ge__(self, other): def __le__(self, other): return self.compare(operator.le, other) + def extract_extrinsic_orientation(self, o): + """Extract extrinsic orientation. + + Parameters + ---------- + o : Orientation + Total orientation. + + Returns + ------- + Orientation + Extrinsic orientation. + + Notes + ----- + The difinition of orientations used here must be consistent with + that used in make_entity_permutations_tensorproduct. + + """ + if not isinstance(o, Orientation): + raise TypeError(f"Expecting an instance of Orientation : got {o}") + dim = len(self.cells) + size_io = 2 # Number of possible intrinsic orientations along each axis. + return o // size_io**dim + + def extract_intrinsic_orientation(self, o, axis): + """Extract intrinsic orientation. + + Parameters + ---------- + o : Orientation + Total orientation. ``//`` and ``%`` must be overloaded in type(o). + axis : int + Reference cell axis for which intrinsic orientation is computed. + + Returns + ------- + Orientation + Intrinsic orientation. + + Notes + ----- + Must be consistent with make_entity_permutations_tensorproduct. + + """ + if not isinstance(o, Orientation): + raise TypeError(f"Expecting an instance of Orientation : got {o}") + dim = len(self.cells) + if axis >= dim: + raise ValueError(f"Must give 0 <= axis < {dim} : got {axis}") + size_io = 2 # Number of possible intrinsic orientations along each axis. + return o % size_io**dim // size_io**(dim - 1 - axis) % size_io + + @property + def extrinsic_orientation_permutation_map(self): + """A map from extrinsic orientations to corresponding axis permutation matrices. + + Notes + ----- + result[eo] gives the physical axis-reference axis permutation matrix corresponding to + eo (extrinsic orientation). + + """ + dim = len(self.cells) + a = numpy.zeros((factorial(dim), dim, dim), dtype=int) + ai = numpy.array(list(make_entity_permutations_simplex(dim - 1, 2).values()), dtype=int).reshape((factorial(dim), dim, 1)) + numpy.put_along_axis(a, ai, 1, axis=2) + return a + class UFCQuadrilateral(Cell): r"""This is the reference quadrilateral with vertices From 98fc7595b5d4ef048b40f04db35b245979fa5b34 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 6 Nov 2024 23:30:29 +0000 Subject: [PATCH 12/13] Fix serendipitous supersmoothness assumption (#102) * Fix serendipitous supersmoothness assumption * unify svd computations --- FIAT/alfeld_sorokina.py | 6 ++---- FIAT/christiansen_hu.py | 6 ++---- FIAT/macro.py | 25 ++++++++++--------------- FIAT/polynomial_set.py | 29 +++++++++++++---------------- 4 files changed, 27 insertions(+), 39 deletions(-) diff --git a/FIAT/alfeld_sorokina.py b/FIAT/alfeld_sorokina.py index 5d3244b78..e09f26d17 100644 --- a/FIAT/alfeld_sorokina.py +++ b/FIAT/alfeld_sorokina.py @@ -41,10 +41,8 @@ def AlfeldSorokinaSpace(ref_el, degree): if len(rows) > 0: dual_mat = numpy.vstack(rows) - _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) - tol = sig[0] * 1E-10 - num_sv = len([s for s in sig if abs(s) > tol]) - coeffs = numpy.tensordot(vt[num_sv:], coeffs, axes=(-1, 0)) + nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) + coeffs = numpy.tensordot(nsp, coeffs, axes=(-1, 0)) return polynomial_set.PolynomialSet(ref_complex, degree, degree, expansion_set, coeffs) diff --git a/FIAT/christiansen_hu.py b/FIAT/christiansen_hu.py index 476ea914d..64fe89659 100644 --- a/FIAT/christiansen_hu.py +++ b/FIAT/christiansen_hu.py @@ -28,10 +28,8 @@ def ChristiansenHuSpace(ref_el, degree, reduced=False): tab = C0.tabulate(Q.get_points(), 1) divC0 = sum(tab[alpha][:, alpha.index(1), :] for alpha in tab if sum(alpha) == 1) - _, sig, vt = numpy.linalg.svd(divC0.T, full_matrices=True) - tol = sig[0] * 1E-10 - num_sv = len([s for s in sig if abs(s) > tol]) - coeffs = numpy.tensordot(vt[num_sv:], C0.get_coeffs(), axes=(-1, 0)) + nsp = polynomial_set.spanning_basis(divC0.T, nullspace=True) + coeffs = numpy.tensordot(nsp, C0.get_coeffs(), axes=(-1, 0)) verts = numpy.array(ref_complex.get_vertices()) WT = verts[-1] diff --git a/FIAT/macro.py b/FIAT/macro.py index f48298c70..77ab89382 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -429,25 +429,21 @@ def __init__(self, ref_el, degree, order=1, vorder=0, shape=(), **kwargs): if len(rows) > 0: dual_mat = numpy.vstack(rows) - _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) - tol = sig[0] * 1E-10 - num_sv = len([s for s in sig if abs(s) > tol]) - coeffs = vt[num_sv:] + coeffs = polynomial_set.spanning_basis(dual_mat, nullspace=True) else: coeffs = numpy.eye(expansion_set.get_num_members(degree)) - if vorder > order + 1: - # Impose C^vorder super-smoothness at interior vertices - # C^order automatically gives C^{order+1} at the interior vertex + # Impose C^vorder super-smoothness at interior vertices + # C^order automatically gives C^{order+dim-1} at the interior vertex + sorder = order + sd - 1 + if vorder > sorder: verts = ref_el.get_vertices() points = [verts[i] for i in ref_el.get_interior_facets(0)] jumps = expansion_set.tabulate_jumps(degree, points, order=vorder) - for r in range(order+2, vorder+1): + for r in range(sorder+1, vorder+1): dual_mat = numpy.dot(numpy.vstack(jumps[r].T), coeffs.T) - _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) - tol = sig[0] * 1E-10 - num_sv = len([s for s in sig if abs(s) > tol]) - coeffs = numpy.dot(vt[num_sv:], coeffs) + nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) + coeffs = numpy.dot(nsp, coeffs) if shape != (): m, n = coeffs.shape @@ -498,8 +494,7 @@ def __init__(self, ref_el, degree, order=0, **kwargs): if len(rows) > 0: dual_mat = numpy.vstack(rows) - _, sig, vt = numpy.linalg.svd(dual_mat, full_matrices=True) - num_sv = len([s for s in sig if abs(s) > 1.e-10]) - coeffs = numpy.tensordot(vt[num_sv:], coeffs, axes=(1, 0)) + nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) + coeffs = numpy.tensordot(nsp, coeffs, axes=(1, 0)) super().__init__(ref_el, degree, degree, expansion_set, coeffs) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 579ec3614..3a3e9ce3a 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -162,28 +162,25 @@ def form_matrix_product(mats, alpha): return result +def spanning_basis(A, nullspace=False, rtol=1e-10): + """Construct a basis that spans the rows of A via SVD. + """ + Aflat = A.reshape(A.shape[0], -1) + u, sig, vt = numpy.linalg.svd(Aflat, full_matrices=True) + atol = rtol * (sig[0] + 1) + num_sv = len([s for s in sig if abs(s) > atol]) + basis = vt[num_sv:] if nullspace else vt[:num_sv] + return numpy.reshape(basis, (-1, *A.shape[1:])) + + def polynomial_set_union_normalized(A, B): """Given polynomial sets A and B, constructs a new polynomial set whose span is the same as that of span(A) union span(B). It may not contain any of the same members of the set, as we construct a span via SVD. """ - new_coeffs = numpy.array(list(A.coeffs) + list(B.coeffs)) - func_shape = new_coeffs.shape[1:] - if len(func_shape) == 1: - (u, sig, vt) = numpy.linalg.svd(new_coeffs) - num_sv = len([s for s in sig if abs(s) > 1.e-10]) - coeffs = vt[:num_sv] - else: - new_shape0 = new_coeffs.shape[0] - new_shape1 = numpy.prod(func_shape) - newshape = (new_shape0, new_shape1) - nc = numpy.reshape(new_coeffs, newshape) - (u, sig, vt) = numpy.linalg.svd(nc, 1) - num_sv = len([s for s in sig if abs(s) > 1.e-10]) - - coeffs = numpy.reshape(vt[:num_sv], (num_sv,) + func_shape) - + new_coeffs = numpy.concatenate((A.coeffs, B.coeffs), axis=0) + coeffs = spanning_basis(new_coeffs) return PolynomialSet(A.get_reference_element(), A.get_degree(), A.get_embedded_degree(), From d167922bf1741e468c74d904eee885cc54f42da4 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 6 Nov 2024 23:30:58 +0000 Subject: [PATCH 13/13] Add super-entities to cell connectivity (#103) --- FIAT/reference_element.py | 59 ++++++++++++++++++++++----------------- 1 file changed, 33 insertions(+), 26 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 95f0d8f8a..822f48de1 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -129,7 +129,7 @@ def linalg_subspace_intersection(A, B): class Cell(object): """Abstract class for a reference cell. Provides accessors for geometry (vertex coordinates) as well as topology (orderings of - vertices that make up edges, facecs, etc.""" + vertices that make up edges, faces, etc.""" def __init__(self, shape, vertices, topology): """The constructor takes a shape code, the physical vertices expressed as a list of tuples of numbers, and the topology of a cell. @@ -160,20 +160,26 @@ def __init__(self, shape, vertices, topology): # Sort for the sake of determinism and by UFC conventions self.sub_entities[dim][e] = sorted(sub_entities) + # Build super-entity dictionary by inverting the sub-entity dictionary + self.super_entities = {dim: {entity: [] for entity in topology[dim]} for dim in topology} + for dim0 in topology: + for e0 in topology[dim0]: + for dim1, e1 in self.sub_entities[dim0][e0]: + self.super_entities[dim1][e1].append((dim0, e0)) + # Build connectivity dictionary for easier queries self.connectivity = {} - for dim0, sub_entities in self.sub_entities.items(): - - # Skip tensor product entities - # TODO: Can we do something better? - if isinstance(dim0, tuple): - continue - - for entity, sub_sub_entities in sorted(sub_entities.items()): - for dim1 in range(dim0+1): - d01_entities = filter(lambda x: x[0] == dim1, sub_sub_entities) - d01_entities = tuple(x[1] for x in d01_entities) - self.connectivity.setdefault((dim0, dim1), []).append(d01_entities) + for dim0 in sorted(topology): + for dim1 in sorted(topology): + self.connectivity[(dim0, dim1)] = [] + + for entity in sorted(topology[dim0]): + children = self.sub_entities[dim0][entity] + parents = self.super_entities[dim0][entity] + for dim1 in sorted(topology): + neighbors = children if dim1 < dim0 else parents + d01_entities = tuple(e for d, e in neighbors if d == dim1) + self.connectivity[(dim0, dim1)].append(d01_entities) def _key(self): """Hashable object key data (excluding type).""" @@ -589,26 +595,27 @@ def get_dimension(self): return self.get_spatial_dimension() def compute_barycentric_coordinates(self, points, entity=None, rescale=False): - """Returns the barycentric coordinates of a list of points on an - entity.""" + """Returns the barycentric coordinates of a list of points on the complex.""" if len(points) == 0: return points if entity is None: entity = (self.get_spatial_dimension(), 0) entity_dim, entity_id = entity top = self.get_topology() - verts = self.get_vertices_of_subcomplex(top[entity_dim][entity_id]) - if entity_dim == self.get_spatial_dimension(): - ref_verts = numpy.eye(entity_dim + 1) - A, b = make_affine_mapping(verts, ref_verts) - else: - v = numpy.transpose(verts) - v = v[:, 1:] - v[:, :1] - A = numpy.linalg.solve(numpy.dot(v.T, v), v.T) - b = -numpy.dot(A, verts[0]) - A = numpy.vstack((-numpy.sum(A, axis=0), A)) - b = numpy.hstack((1 - numpy.sum(b, axis=0), b)) + 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] + + cell_verts = self.get_vertices_of_subcomplex(subcomplex) + ref_verts = numpy.eye(sd + 1) + 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 h = 1 / numpy.linalg.norm(A, axis=1)