diff --git a/FIAT/__init__.py b/FIAT/__init__.py index b657b893f..cc399210a 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -7,6 +7,7 @@ # Import finite element classes from FIAT.finite_element import FiniteElement, CiarletElement # noqa: F401 from FIAT.argyris import Argyris +from FIAT.bernstein import Bernstein from FIAT.bell import Bell from FIAT.argyris import QuinticArgyris from FIAT.brezzi_douglas_marini import BrezziDouglasMarini @@ -14,6 +15,7 @@ from FIAT.discontinuous_lagrange import DiscontinuousLagrange from FIAT.discontinuous_taylor import DiscontinuousTaylor from FIAT.discontinuous_raviart_thomas import DiscontinuousRaviartThomas +from FIAT.discontinuous_pc import DPC from FIAT.hermite import CubicHermite from FIAT.lagrange import Lagrange from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre @@ -47,12 +49,14 @@ # List of supported elements and mapping to element classes supported_elements = {"Argyris": Argyris, "Bell": Bell, + "Bernstein": Bernstein, "Brezzi-Douglas-Marini": BrezziDouglasMarini, "Brezzi-Douglas-Fortin-Marini": BrezziDouglasFortinMarini, "Bubble": Bubble, "FacetBubble": FacetBubble, "Crouzeix-Raviart": CrouzeixRaviart, "Discontinuous Lagrange": DiscontinuousLagrange, + "DPC": DPC, "Discontinuous Taylor": DiscontinuousTaylor, "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas, "Hermite": CubicHermite, diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py new file mode 100644 index 000000000..251f94204 --- /dev/null +++ b/FIAT/bernstein.py @@ -0,0 +1,199 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2018 Miklós Homolya +# +# This file is part of FIAT. +# +# FIAT is free software: you can redistribute it and/or modify it +# under the terms of the GNU Lesser General Public License as +# published by the Free Software Foundation, either version 3 of the +# License, or (at your option) any later version. +# +# FIAT is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +# License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with FIAT. If not, see . + +import math +import numpy + +from FIAT.finite_element import FiniteElement +from FIAT.dual_set import DualSet +from FIAT.polynomial_set import mis + + +class BernsteinDualSet(DualSet): + """The dual basis for Bernstein elements.""" + + def __init__(self, ref_el, degree): + # Initialise data structures + topology = ref_el.get_topology() + entity_ids = {dim: {entity_i: [] + for entity_i in entities} + for dim, entities in topology.items()} + + # Calculate inverse topology + inverse_topology = {vertices: (dim, entity_i) + for dim, entities in topology.items() + for entity_i, vertices in entities.items()} + + # Generate triangular barycentric indices + dim = ref_el.get_spatial_dimension() + kss = mis(dim + 1, degree) + + # Fill data structures + nodes = [] + for i, ks in enumerate(kss): + vertices, = numpy.nonzero(ks) + entity_dim, entity_i = inverse_topology[tuple(vertices)] + entity_ids[entity_dim][entity_i].append(i) + + # Leave nodes unimplemented for now + nodes.append(None) + + super(BernsteinDualSet, self).__init__(nodes, ref_el, entity_ids) + + +class Bernstein(FiniteElement): + """A finite element with Bernstein polynomials as basis functions.""" + + def __init__(self, ref_el, degree): + dual = BernsteinDualSet(ref_el, degree) + k = 0 # 0-form + super(Bernstein, self).__init__(ref_el, dual, degree, k) + + def degree(self): + """The degree of the polynomial space.""" + return self.get_order() + + def value_shape(self): + """The value shape of the finite element functions.""" + return () + + def tabulate(self, order, points, entity=None): + """Return tabulated values of derivatives up to given order of + basis functions at given points. + + :arg order: The maximum order of derivative. + :arg points: An iterable of points. + :arg entity: Optional (dimension, entity number) pair + indicating which topological entity of the + reference element to tabulate on. If ``None``, + default cell-wise tabulation is performed. + """ + # Transform points to reference cell coordinates + ref_el = self.get_reference_element() + if entity is None: + entity = (ref_el.get_spatial_dimension(), 0) + + entity_dim, entity_id = entity + entity_transform = ref_el.get_entity_transform(entity_dim, entity_id) + cell_points = list(map(entity_transform, points)) + + # Construct Cartesian to Barycentric coordinate mapping + vs = numpy.asarray(ref_el.get_vertices()) + B2R = numpy.vstack([vs.T, numpy.ones(len(vs))]) + R2B = numpy.linalg.inv(B2R) + + B = numpy.hstack([cell_points, + numpy.ones((len(cell_points), 1))]).dot(R2B.T) + + # Evaluate everything + deg = self.degree() + dim = ref_el.get_spatial_dimension() + raw_result = {(alpha, i): vec + for i, ks in enumerate(mis(dim + 1, deg)) + for o in range(order + 1) + for alpha, vec in bernstein_Dx(B, ks, o, R2B).items()} + + # Rearrange result + space_dim = self.space_dimension() + dtype = numpy.array(list(raw_result.values())).dtype + result = {alpha: numpy.zeros((space_dim, len(cell_points)), dtype=dtype) + for o in range(order + 1) + for alpha in mis(dim, o)} + for (alpha, i), vec in raw_result.items(): + result[alpha][i, :] = vec + return result + + +def bernstein_db(points, ks, alpha=None): + """Evaluates Bernstein polynomials or its derivative at barycentric + points. + + :arg points: array of points in barycentric coordinates + :arg ks: exponents defining the Bernstein polynomial + :arg alpha: derivative tuple + + :returns: array of Bernstein polynomial values at given points. + """ + points = numpy.asarray(points) + ks = numpy.array(tuple(ks)) + + N, d_1 = points.shape + assert d_1 == len(ks) + + if alpha is None: + alpha = numpy.zeros(d_1) + else: + alpha = numpy.array(tuple(alpha)) + assert d_1 == len(alpha) + + ls = ks - alpha + if any(k < 0 for k in ls): + return numpy.zeros(len(points)) + elif all(k == 0 for k in ls): + return numpy.ones(len(points)) + else: + # Calculate coefficient + coeff = math.factorial(ks.sum()) + for k in ls: + coeff //= math.factorial(k) + return coeff * numpy.prod(points**ls, axis=1) + + +def bernstein_Dx(points, ks, order, R2B): + """Evaluates Bernstein polynomials or its derivatives according to + reference coordinates. + + :arg points: array of points in BARYCENTRIC COORDINATES + :arg ks: exponents defining the Bernstein polynomial + :arg alpha: derivative order (returns all derivatives of this + specified order) + :arg R2B: linear mapping from reference to barycentric coordinates + + :returns: dictionary mapping from derivative tuples to arrays of + Bernstein polynomial values at given points. + """ + points = numpy.asarray(points) + ks = tuple(ks) + + N, d_1 = points.shape + assert d_1 == len(ks) + + # Collect derivatives according to barycentric coordinates + Db_map = {alpha: bernstein_db(points, ks, alpha) + for alpha in mis(d_1, order)} + + # Arrange derivative tensor (barycentric coordinates) + dtype = numpy.array(list(Db_map.values())).dtype + Db_shape = (d_1,) * order + Db_tensor = numpy.empty(Db_shape + (N,), dtype=dtype) + for ds in numpy.ndindex(Db_shape): + alpha = [0] * d_1 + for d in ds: + alpha[d] += 1 + Db_tensor[ds + (slice(None),)] = Db_map[tuple(alpha)] + + # Coordinate transformation: barycentric -> reference + result = {} + for alpha in mis(d_1 - 1, order): + values = Db_tensor + for d, k in enumerate(alpha): + for _ in range(k): + values = R2B[:, d].dot(values) + result[alpha] = values + return result diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py new file mode 100644 index 000000000..8c5683109 --- /dev/null +++ b/FIAT/discontinuous_pc.py @@ -0,0 +1,119 @@ +# Copyright (C) 2018 Cyrus Cheng (Imperial College London) +# +# This file is part of FIAT. +# +# FIAT is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# FIAT is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with FIAT. If not, see . +# +# Modified by David A. Ham (david.ham@imperial.ac.uk), 2018 + +from FIAT import finite_element, polynomial_set, dual_set, functional +from FIAT.P0 import P0Dual +from FIAT.reference_element import (Point, + DefaultLine, + UFCInterval, + UFCQuadrilateral, + UFCHexahedron, + UFCTriangle, + UFCTetrahedron, + make_affine_mapping) +import numpy as np + +hypercube_simplex_map = {Point(): Point(), + DefaultLine(): DefaultLine(), + UFCInterval(): UFCInterval(), + UFCQuadrilateral(): UFCTriangle(), + UFCHexahedron(): UFCTetrahedron()} + + +class DPC0(finite_element.CiarletElement): + def __init__(self, ref_el): + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) + dual = P0Dual(ref_el) + degree = 0 + formdegree = ref_el.get_spatial_dimension() # n-form + super(DPC0, self).__init__(poly_set=poly_set, + dual=dual, + order=degree, + ref_el=ref_el, + formdegree=formdegree) + + +class DPCDualSet(dual_set.DualSet): + """The dual basis for DPC elements. This class works for + hypercubes of any dimension. Nodes are point evaluation at + equispaced points. This is the discontinuous version where + all nodes are topologically associated with the cell itself""" + + def __init__(self, ref_el, degree): + entity_ids = {} + nodes = [] + + # Change coordinates here. + # Vertices of the simplex corresponding to the reference element. + v_simplex = hypercube_simplex_map[ref_el].get_vertices() + # Vertices of the reference element. + v_hypercube = ref_el.get_vertices() + # For the mapping, first two vertices are unchanged in all dimensions. + v_ = [v_hypercube[0], v_hypercube[int(-0.5*len(v_hypercube))]] + + # For dimension 1 upwards, + # take the next vertex and map it to the midpoint of the edge/face it belongs to, and shares + # with no other points. + for d in range(1, ref_el.get_dimension()): + v_.append(tuple(np.asarray(v_hypercube[ref_el.get_dimension() - d] + + np.average(np.asarray(v_hypercube[::2]), axis=0)))) + A, b = make_affine_mapping(v_simplex, tuple(v_)) # Make affine mapping to be used later. + + # make nodes by getting points + # need to do this dimension-by-dimension, facet-by-facet + top = hypercube_simplex_map[ref_el].get_topology() + cube_topology = ref_el.get_topology() + + cur = 0 + for dim in sorted(top): + entity_ids[dim] = {} + for entity in sorted(top[dim]): + pts_cur = hypercube_simplex_map[ref_el].make_points(dim, entity, degree) + pts_cur = [tuple(np.matmul(A, np.array(x)) + b) for x in pts_cur] + nodes_cur = [functional.PointEvaluation(ref_el, x) + for x in pts_cur] + nnodes_cur = len(nodes_cur) + nodes += nodes_cur + cur += nnodes_cur + for entity in sorted(cube_topology[dim]): + entity_ids[dim][entity] = [] + + entity_ids[dim][0] = list(range(len(nodes))) + super(DPCDualSet, self).__init__(nodes, ref_el, entity_ids) + + +class HigherOrderDPC(finite_element.CiarletElement): + """The DPC finite element. It is what it is.""" + + def __init__(self, ref_el, degree): + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], degree) + dual = DPCDualSet(ref_el, degree) + formdegree = ref_el.get_spatial_dimension() # n-form + super(HigherOrderDPC, self).__init__(poly_set=poly_set, + dual=dual, + order=degree, + ref_el=ref_el, + formdegree=formdegree) + + +def DPC(ref_el, degree): + if degree == 0: + return DPC0(ref_el) + else: + return HigherOrderDPC(ref_el, degree) diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index 7ca937534..15ead7129 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -120,8 +120,8 @@ class CiarletElement(FiniteElement): basis generated from polynomials encoded in a `PolynomialSet`. """ - def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine"): - ref_el = poly_set.get_reference_element() + def __init__(self, poly_set, dual, order, formdegree=None, mapping="affine", ref_el=None): + ref_el = ref_el or poly_set.get_reference_element() super(CiarletElement, self).__init__(ref_el, dual, order, formdegree, mapping) # build generalized Vandermonde matrix diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index dde620a2d..3b9353c21 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -348,19 +348,13 @@ def barycentric_coordinates(points, vertices): """ # Form mapping matrix - last = np.asarray(vertices[-1]) - T = np.matrix([np.array(v) - last for v in vertices[:-1]]).T + T = (np.asarray(vertices[:-1]) - vertices[-1]).T invT = np.linalg.inv(T) - # Compute the barycentric coordinates for all points - coords = [] - for p in points: - y = np.asarray(p) - last - bary = invT.dot(y.T) - bary = [bary[(0, i)] for i in range(len(y))] - bary += [1.0 - sum(bary)] - coords.append(bary) - return coords + points = np.asarray(points) + bary = np.einsum("ij,kj->ki", invT, (points - vertices[-1])) + last = (1 - bary.sum(axis=1)) + return np.concatenate([bary, last[..., np.newaxis]], axis=1) def map_from_reference_facet(point, vertices): diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 3198a028a..fcd264df6 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -29,7 +29,7 @@ class NedelecSecondKindDual(DualSet): - """ + r""" This class represents the dual basis for the Nedelec H(curl) elements of the second kind. The degrees of freedom (L) for the elements of the k'th degree are @@ -160,11 +160,11 @@ def _generate_face_dofs(self, cell, degree, offset): # Map Phis -> phis (reference values to physical values) J = Q_face.jacobian() - scale = 1.0 / numpy.sqrt(numpy.linalg.det(J.transpose() * J)) + scale = 1.0 / numpy.sqrt(numpy.linalg.det(numpy.dot(J.T, J))) phis = numpy.ndarray((d, num_quad_points)) for i in range(num_rts): for q in range(num_quad_points): - phi_i_q = scale * J * numpy.matrix(Phis[i, :, q]).transpose() + phi_i_q = scale * numpy.dot(J, Phis[numpy.newaxis, i, :, q].T) for j in range(d): phis[j, q] = phi_i_q[j] diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py index 7b8b5490c..95292bb0b 100644 --- a/FIAT/quadrature.py +++ b/FIAT/quadrature.py @@ -211,15 +211,12 @@ def __init__(self, face_number, degree): # Use tet to map points and weights on the appropriate face vertices = [numpy.array(list(vertex)) for vertex in vertices] x0 = vertices[0] - J = numpy.matrix([vertices[1] - x0, vertices[2] - x0]).transpose() - x0 = numpy.matrix(x0).transpose() + J = numpy.vstack([vertices[1] - x0, vertices[2] - x0]).T # This is just a very numpyfied way of writing J*p + x0: - F = lambda p: \ - numpy.array(J*numpy.matrix(p).transpose() + x0).flatten() - points = numpy.array([F(p) for p in ref_points]) + points = numpy.einsum("ij,kj->ki", J, ref_points) + x0 # Map weights: multiply reference weights by sqrt(|J^T J|) - detJTJ = numpy.linalg.det(J.transpose() * J) + detJTJ = numpy.linalg.det(numpy.dot(J.T, J)) weights = numpy.sqrt(detJTJ) * ref_weights # Initialize super class with new points and weights @@ -255,6 +252,14 @@ def make_quadrature(ref_el, m): return CollapsedQuadratureTriangleRule(ref_el, m) elif ref_el.get_shape() == reference_element.TETRAHEDRON: return CollapsedQuadratureTetrahedronRule(ref_el, m) + elif ref_el.get_shape() == reference_element.QUADRILATERAL: + line_rule = GaussJacobiQuadratureLineRule(ref_el.construct_subelement(1), m) + return make_tensor_product_quadrature(line_rule, line_rule) + elif ref_el.get_shape() == reference_element.HEXAHEDRON: + line_rule = GaussJacobiQuadratureLineRule(ref_el.construct_subelement(1), m) + return make_tensor_product_quadrature(line_rule, line_rule, line_rule) + else: + raise ValueError("Unable to make quadrature for cell: %s" % ref_el) def make_tensor_product_quadrature(*quad_rules): diff --git a/setup.cfg b/setup.cfg index bcaf19f1d..82aba20a5 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,4 +1,4 @@ [flake8] -ignore = E501,E226,E731 +ignore = E501,E226,E731,W504 exclude = .git,__pycache__,doc/sphinx/source/conf.py,build,dist min-version = 3.0 diff --git a/test/unit/test_bernstein.py b/test/unit/test_bernstein.py new file mode 100644 index 000000000..2b4714a8a --- /dev/null +++ b/test/unit/test_bernstein.py @@ -0,0 +1,85 @@ +# -*- coding: utf-8 -*- +# +# Copyright (C) 2018 Miklós Homolya +# +# This file is part of FIAT. +# +# FIAT is free software: you can redistribute it and/or modify it +# under the terms of the GNU Lesser General Public License as +# published by the Free Software Foundation, either version 3 of the +# License, or (at your option) any later version. +# +# FIAT is distributed in the hope that it will be useful, but WITHOUT +# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +# License for more details. +# +# You should have received a copy of the GNU Lesser General Public +# License along with FIAT. If not, see . + +import numpy +import pytest + +from FIAT.reference_element import ufc_simplex +from FIAT.bernstein import Bernstein +from FIAT.quadrature_schemes import create_quadrature + + +D02 = numpy.array([ + [0.65423405, 1.39160021, 0.65423405, 3.95416573, 1.39160021, 3.95416573], + [3.95416573, 3.95416573, 1.39160021, 1.39160021, 0.65423405, 0.65423405], + [0.0831321, -2.12896637, 2.64569763, -7.25409741, 1.17096531, -6.51673126], + [0., 0., 0., 0., 0., 0.], + [-7.90833147, -7.90833147, -2.78320042, -2.78320042, -1.30846811, -1.30846811], + [-2.12896637, 0.0831321, -7.25409741, 2.64569763, -6.51673126, 1.17096531], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.], + [3.95416573, 3.95416573, 1.39160021, 1.39160021, 0.65423405, 0.65423405], + [1.39160021, 0.65423405, 3.95416573, 0.65423405, 3.95416573, 1.39160021], +]) + +D11 = numpy.array([ + [0.65423405, 1.39160021, 0.65423405, 3.95416573, 1.39160021, 3.95416573], + [3.29993168, 2.56256552, 0.73736616, -2.56256552, -0.73736616, -3.29993168], + [0.73736616, -0.73736616, 3.29993168, -3.29993168, 2.56256552, -2.56256552], + [-3.95416573, -3.95416573, -1.39160021, -1.39160021, -0.65423405, -0.65423405], + [-4.69153189, -3.21679958, -4.69153189, 1.90833147, -3.21679958, 1.90833147], + [-1.39160021, -0.65423405, -3.95416573, -0.65423405, -3.95416573, -1.39160021], + [0., 0., 0., 0., 0., 0.], + [3.95416573, 3.95416573, 1.39160021, 1.39160021, 0.65423405, 0.65423405], + [1.39160021, 0.65423405, 3.95416573, 0.65423405, 3.95416573, 1.39160021], + [0., 0., 0., 0., 0., 0.], +]) + +D20 = numpy.array([ + [0.65423405, 1.39160021, 0.65423405, 3.95416573, 1.39160021, 3.95416573], + [2.64569763, 1.17096531, 0.0831321, -6.51673126, -2.12896637, -7.25409741], + [1.39160021, 0.65423405, 3.95416573, 0.65423405, 3.95416573, 1.39160021], + [-7.25409741, -6.51673126, -2.12896637, 1.17096531, 0.0831321, 2.64569763], + [-2.78320042, -1.30846811, -7.90833147, -1.30846811, -7.90833147, -2.78320042], + [0., 0., 0., 0., 0., 0.], + [3.95416573, 3.95416573, 1.39160021, 1.39160021, 0.65423405, 0.65423405], + [1.39160021, 0.65423405, 3.95416573, 0.65423405, 3.95416573, 1.39160021], + [0., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.], +]) + + +def test_bernstein_2nd_derivatives(): + ref_el = ufc_simplex(2) + degree = 3 + + elem = Bernstein(ref_el, degree) + rule = create_quadrature(ref_el, degree) + points = rule.get_points() + + actual = elem.tabulate(2, points) + + assert numpy.allclose(D02, actual[(0, 2)]) + assert numpy.allclose(D11, actual[(1, 1)]) + assert numpy.allclose(D20, actual[(2, 0)]) + + +if __name__ == '__main__': + import os + pytest.main(os.path.abspath(__file__)) diff --git a/test/unit/test_discontinuous_pc.py b/test/unit/test_discontinuous_pc.py new file mode 100644 index 000000000..afb29b065 --- /dev/null +++ b/test/unit/test_discontinuous_pc.py @@ -0,0 +1,51 @@ +# Copyright (C) 2016 Imperial College London and others +# +# This file is part of FIAT. +# +# FIAT is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# FIAT is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with FIAT. If not, see . +# +# Authors: +# +# David Ham + +import pytest +import numpy as np + + +@pytest.mark.parametrize("dim, degree", [(dim, degree) + for dim in range(1, 4) + for degree in range(4)]) +def test_basis_values(dim, degree): + """Ensure that integrating a simple monomial produces the expected results.""" + from FIAT import ufc_cell, make_quadrature + from FIAT.discontinuous_pc import DPC + + cell = np.array([None, 'interval', 'quadrilateral', 'hexahedron']) + s = ufc_cell(cell[dim]) + q = make_quadrature(s, degree + 1) + + fe = DPC(s, degree) + tab = fe.tabulate(0, q.pts)[(0,) * dim] + + for test_degree in range(degree + 1): + coefs = [n(lambda x: x[0]**test_degree) for n in fe.dual.nodes] + integral = np.float(np.dot(coefs, np.dot(tab, q.wts))) + reference = np.dot([x[0]**test_degree + for x in q.pts], q.wts) + assert np.isclose(integral, reference, rtol=1e-14) + + +if __name__ == '__main__': + import os + pytest.main(os.path.abspath(__file__))