From 812395ae3f7daa1c2b1b374ddf4d0794ba07393b Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Wed, 3 Oct 2018 22:54:35 +0100 Subject: [PATCH 01/46] Remove reliance on numpy.matrix Squashes more deprecation warnings. --- FIAT/hdiv_trace.py | 16 +++++----------- FIAT/nedelec_second_kind.py | 4 ++-- FIAT/quadrature.py | 9 +++------ 3 files changed, 10 insertions(+), 19 deletions(-) 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..1c53842c0 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -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..035326093 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 From 691d3fa7c09a0afeb88afca772a4b2fa4d8e9f38 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Thu, 25 Oct 2018 09:05:22 +0200 Subject: [PATCH 02/46] fix new flake8 rules --- FIAT/nedelec_second_kind.py | 2 +- setup.cfg | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 3198a028a..81a34ded3 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 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 From 77487342093f8039256822d9f2b53b885e8e77f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Wed, 10 Oct 2018 21:48:41 +0200 Subject: [PATCH 03/46] add incomplete Bernstein element No derivative tabulations yet. --- FIAT/__init__.py | 2 + FIAT/bernstein.py | 117 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 119 insertions(+) create mode 100644 FIAT/bernstein.py diff --git a/FIAT/__init__.py b/FIAT/__init__.py index b657b893f..8e60bb270 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 @@ -47,6 +48,7 @@ # 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, diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py new file mode 100644 index 000000000..c42ab7b2e --- /dev/null +++ b/FIAT/bernstein.py @@ -0,0 +1,117 @@ +# -*- 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 General Public License +# along with FIAT. If not, see . + +import itertools +import numpy + +from FIAT.finite_element import FiniteElement +from FIAT.dual_set import DualSet + + +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() + alphas = [(degree - sum(beta),) + tuple(reversed(beta)) + for beta in itertools.product(range(degree + 1), repeat=dim) + if sum(beta) <= degree] + + # Fill data structures + nodes = [] + for i, alpha in enumerate(alphas): + vertices, = numpy.nonzero(alpha) + 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. + """ + # Retrieve entity transformation + 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) + + # Construct Cartesian to Barycentric coordinate mapping + B2C = numpy.hstack([ref_el.get_vertices(), + numpy.ones((len(ref_el.get_vertices()), 1))]) + C2B = numpy.linalg.inv(B2C) + + # Array of barycentric point coordinates + Bs = numpy.array([tuple(entity_transform(point)) + (1,) + for point in points]).dot(C2B) + + # Generate triangular barycentric indices + dim = ref_el.get_spatial_dimension() + deg = self.degree() + alphas = [(deg - sum(beta),) + tuple(reversed(beta)) + for beta in itertools.product(range(deg + 1), repeat=dim) + if sum(beta) <= deg] + + assert order == 0 + result = numpy.zeros((len(alphas), len(Bs))) + for i, alpha in enumerate(alphas): + for j, bs in enumerate(Bs): + result[i, j] = numpy.prod([b**k for b, k in zip(bs, alpha)]) + return {(0,) * ref_el.get_spatial_dimension(): result} From cfd0e4336d9384f96f034c678f8638ea13761228 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Thu, 11 Oct 2018 10:35:46 +0200 Subject: [PATCH 04/46] add forgotten scaling factor --- FIAT/bernstein.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index c42ab7b2e..06476c807 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -18,6 +18,7 @@ # along with FIAT. If not, see . import itertools +import math import numpy from FIAT.finite_element import FiniteElement @@ -113,5 +114,8 @@ def tabulate(self, order, points, entity=None): result = numpy.zeros((len(alphas), len(Bs))) for i, alpha in enumerate(alphas): for j, bs in enumerate(Bs): - result[i, j] = numpy.prod([b**k for b, k in zip(bs, alpha)]) + c = math.factorial(deg) + for k in alpha: + c = c // math.factorial(k) + result[i, j] = c * numpy.prod([b**k for b, k in zip(bs, alpha)]) return {(0,) * ref_el.get_spatial_dimension(): result} From dd808d34ec1397e0cdddb8779f5b1b7f60ab5e86 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Mon, 15 Oct 2018 08:46:42 +0200 Subject: [PATCH 05/46] add very slow implementation of derivatives --- FIAT/bernstein.py | 50 +++++++++++++++++++++++++++-------------------- 1 file changed, 29 insertions(+), 21 deletions(-) diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index 06476c807..476ab5c82 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -19,10 +19,13 @@ import itertools import math + import numpy +import sympy from FIAT.finite_element import FiniteElement from FIAT.dual_set import DualSet +from FIAT.polynomial_set import mis class BernsteinDualSet(DualSet): @@ -86,36 +89,41 @@ def tabulate(self, order, points, entity=None): reference element to tabulate on. If ``None``, default cell-wise tabulation is performed. """ - # Retrieve entity transformation + # 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 - B2C = numpy.hstack([ref_el.get_vertices(), - numpy.ones((len(ref_el.get_vertices()), 1))]) - C2B = numpy.linalg.inv(B2C) + vs = numpy.asarray(ref_el.get_vertices()) + B2R = numpy.vstack([vs.T, numpy.ones(len(vs))]) + R2B = numpy.linalg.inv(B2R) - # Array of barycentric point coordinates - Bs = numpy.array([tuple(entity_transform(point)) + (1,) - for point in points]).dot(C2B) + dim = ref_el.get_spatial_dimension() + X = sympy.symbols('X Y Z')[:dim] + B = R2B.dot(X + (1,)) # Generate triangular barycentric indices - dim = ref_el.get_spatial_dimension() deg = self.degree() - alphas = [(deg - sum(beta),) + tuple(reversed(beta)) - for beta in itertools.product(range(deg + 1), repeat=dim) - if sum(beta) <= deg] - - assert order == 0 - result = numpy.zeros((len(alphas), len(Bs))) - for i, alpha in enumerate(alphas): - for j, bs in enumerate(Bs): - c = math.factorial(deg) - for k in alpha: - c = c // math.factorial(k) - result[i, j] = c * numpy.prod([b**k for b, k in zip(bs, alpha)]) - return {(0,) * ref_el.get_spatial_dimension(): result} + etas = [(deg - sum(beta),) + tuple(reversed(beta)) + for beta in itertools.product(range(deg + 1), repeat=dim) + if sum(beta) <= deg] + + result = {} + for D in range(order + 1): + for alpha in mis(dim, D): + table = numpy.zeros((len(etas), len(cell_points))) + for i, eta in enumerate(etas): + for j, point in enumerate(cell_points): + c = math.factorial(deg) + for k in eta: + c = c // math.factorial(k) + b = c * (B ** eta).prod() + e = sympy.diff(b, *zip(X, alpha)) + table[i, j] = e.subs(dict(zip(X, point))).evalf() + result[alpha] = table + return result From 055d3063970a3e1d7ac0c2eccb1d9d7137e0051f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Mon, 22 Oct 2018 11:29:52 +0200 Subject: [PATCH 06/46] test higher derivatives of Bernstein elements --- test/unit/test_bernstein.py | 85 +++++++++++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) create mode 100644 test/unit/test_bernstein.py diff --git a/test/unit/test_bernstein.py b/test/unit/test_bernstein.py new file mode 100644 index 000000000..74de5d03c --- /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 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., 0., 0., 0., 0., 0.], + [0., 0., 0., 0., 0., 0.], + [0.0831321, -2.12896637, 2.64569763, -7.25409741, 1.17096531, -6.51673126], + [-7.90833147, -7.90833147, -2.78320042, -2.78320042, -1.30846811, -1.30846811], + [0., 0., 0., 0., 0., 0.], + [-2.12896637, 0.0831321, -7.25409741, 2.64569763, -6.51673126, 1.17096531], + [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], + [-3.95416573, -3.95416573, -1.39160021, -1.39160021, -0.65423405, -0.65423405], + [0., 0., 0., 0., 0., 0.], + [0.73736616, -0.73736616, 3.29993168, -3.29993168, 2.56256552, -2.56256552], + [-4.69153189, -3.21679958, -4.69153189, 1.90833147, -3.21679958, 1.90833147], + [3.95416573, 3.95416573, 1.39160021, 1.39160021, 0.65423405, 0.65423405], + [-1.39160021, -0.65423405, -3.95416573, -0.65423405, -3.95416573, -1.39160021], + [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], + [-7.25409741, -6.51673126, -2.12896637, 1.17096531, 0.0831321, 2.64569763], + [3.95416573, 3.95416573, 1.39160021, 1.39160021, 0.65423405, 0.65423405], + [1.39160021, 0.65423405, 3.95416573, 0.65423405, 3.95416573, 1.39160021], + [-2.78320042, -1.30846811, -7.90833147, -1.30846811, -7.90833147, -2.78320042], + [1.39160021, 0.65423405, 3.95416573, 0.65423405, 3.95416573, 1.39160021], + [0., 0., 0., 0., 0., 0.], + [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__)) From 0090e88c97c92220d2363d1f8c6dde61ac7a80b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Tue, 23 Oct 2018 09:46:16 +0200 Subject: [PATCH 07/46] speed up Bernstein evaluations Banish expensive sympy and use fast numpy. --- FIAT/bernstein.py | 96 +++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 85 insertions(+), 11 deletions(-) diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index 476ab5c82..0a4a9c9ff 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -21,7 +21,6 @@ import math import numpy -import sympy from FIAT.finite_element import FiniteElement from FIAT.dual_set import DualSet @@ -103,12 +102,14 @@ def tabulate(self, order, points, entity=None): B2R = numpy.vstack([vs.T, numpy.ones(len(vs))]) R2B = numpy.linalg.inv(B2R) - dim = ref_el.get_spatial_dimension() - X = sympy.symbols('X Y Z')[:dim] - B = R2B.dot(X + (1,)) + B = numpy.hstack([cell_points, + numpy.ones((len(cell_points), 1))]).dot(R2B.T) + # X = sympy.symbols('X Y Z')[:dim] + # B = R2B.dot(X + (1,)) # Generate triangular barycentric indices deg = self.degree() + dim = ref_el.get_spatial_dimension() etas = [(deg - sum(beta),) + tuple(reversed(beta)) for beta in itertools.product(range(deg + 1), repeat=dim) if sum(beta) <= deg] @@ -118,12 +119,85 @@ def tabulate(self, order, points, entity=None): for alpha in mis(dim, D): table = numpy.zeros((len(etas), len(cell_points))) for i, eta in enumerate(etas): - for j, point in enumerate(cell_points): - c = math.factorial(deg) - for k in eta: - c = c // math.factorial(k) - b = c * (B ** eta).prod() - e = sympy.diff(b, *zip(X, alpha)) - table[i, j] = e.subs(dict(zip(X, point))).evalf() + table[i, :] = bernstein_dx(B, eta, alpha, R2B) + # for j, point in enumerate(cell_points): + # # c = math.factorial(deg) + # # for k in eta: + # # c = c // math.factorial(k) + # b = c * (B ** eta).prod() + # e = sympy.diff(b, *zip(X, alpha)) + # table[i, j] = e.subs(dict(zip(X, point))).evalf() result[alpha] = table return result + + +def bernstein_b(points, alpha): + """Evaluates Bernstein polynomials at barycentric points. + + :arg points: array of points in barycentric coordinates + :arg alpha: exponents defining the Bernstein polynomial + :returns: array of Bernstein polynomial values at given points. + """ + points = numpy.asarray(points) + alpha = tuple(alpha) + + N, d_1 = points.shape + assert d_1 == len(alpha) + if any(k < 0 for k in alpha): + return numpy.zeros(len(points)) + elif all(k == 0 for k in alpha): + return numpy.ones(len(points)) + else: + c = math.factorial(sum(alpha)) + for k in alpha: + c = c // math.factorial(k) + return c * numpy.prod(points**alpha, axis=1) + + +def bernstein_db(points, alpha, delta): + points = numpy.asarray(points) + alpha = tuple(alpha) + delta = tuple(delta) + + N, d_1 = points.shape + assert d_1 == len(alpha) == len(delta) + + # Calculate derivative factor + c = 1 + for _, i in zip(range(sum(delta)), range(sum(alpha), 0, -1)): + c *= i + + alpha_ = numpy.array(alpha) - numpy.array(delta) + return c * bernstein_b(points, alpha_) + + +def bernstein_Db(points, alpha, order): + points = numpy.asarray(points) + alpha = tuple(alpha) + + N, d_1 = points.shape + assert d_1 == len(alpha) + Dshape = (d_1,) * order + + result = numpy.empty(Dshape + (N,)) + for indices in numpy.ndindex(Dshape): + delta = [0] * d_1 + for i in indices: + delta[i] += 1 + result[indices + (slice(None),)] = bernstein_db(points, alpha, delta) + return result + + +def bernstein_dx(points, alpha, delta, R2B): + points = numpy.asarray(points) + alpha = tuple(alpha) + delta = tuple(delta) + + N, d_1 = points.shape + assert d_1 == len(alpha) == len(delta) + 1 + + result = bernstein_Db(points, alpha, sum(delta)) + for d, c in enumerate(delta): + for _ in range(c): + result = R2B[:, d].dot(result) + return result From 311b42b16fe24817e897dbd7fa45e676def82e89 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Tue, 23 Oct 2018 14:37:09 +0200 Subject: [PATCH 08/46] partially clean up Bernstein evaluations Performance could be improved further, but it is already much faster than Lagrange elements. (Although, Lagrange elements are less relevant for high-order than Bernstein elements.) --- FIAT/bernstein.py | 119 ++++++++++++++++++++-------------------------- 1 file changed, 52 insertions(+), 67 deletions(-) diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index 0a4a9c9ff..39c7dbcc7 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -104,8 +104,6 @@ def tabulate(self, order, points, entity=None): B = numpy.hstack([cell_points, numpy.ones((len(cell_points), 1))]).dot(R2B.T) - # X = sympy.symbols('X Y Z')[:dim] - # B = R2B.dot(X + (1,)) # Generate triangular barycentric indices deg = self.degree() @@ -117,87 +115,74 @@ def tabulate(self, order, points, entity=None): result = {} for D in range(order + 1): for alpha in mis(dim, D): - table = numpy.zeros((len(etas), len(cell_points))) - for i, eta in enumerate(etas): - table[i, :] = bernstein_dx(B, eta, alpha, R2B) - # for j, point in enumerate(cell_points): - # # c = math.factorial(deg) - # # for k in eta: - # # c = c // math.factorial(k) - # b = c * (B ** eta).prod() - # e = sympy.diff(b, *zip(X, alpha)) - # table[i, j] = e.subs(dict(zip(X, point))).evalf() - result[alpha] = table + result[alpha] = numpy.zeros((len(etas), len(cell_points))) + for i, eta in enumerate(etas): + for alpha, vec in bernstein_Dx(B, eta, D, R2B).items(): + result[alpha][i, :] = vec return result -def bernstein_b(points, alpha): - """Evaluates Bernstein polynomials at barycentric points. +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 alpha: exponents defining the Bernstein polynomial + :arg ks: exponents defining the Bernstein polynomial + :arg alpha: derivative tuple + :returns: array of Bernstein polynomial values at given points. """ points = numpy.asarray(points) - alpha = tuple(alpha) + ks = numpy.array(tuple(ks)) N, d_1 = points.shape - assert d_1 == len(alpha) - if any(k < 0 for k in alpha): + 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 alpha): + elif all(k == 0 for k in ls): return numpy.ones(len(points)) else: - c = math.factorial(sum(alpha)) - for k in alpha: - c = c // math.factorial(k) - return c * numpy.prod(points**alpha, axis=1) - - -def bernstein_db(points, alpha, delta): - points = numpy.asarray(points) - alpha = tuple(alpha) - delta = tuple(delta) - - N, d_1 = points.shape - assert d_1 == len(alpha) == len(delta) + # Calculate coefficient + coeff = math.factorial(ks.sum()) + for k in ls: + coeff //= math.factorial(k) + return coeff * numpy.prod(points**ls, axis=1) - # Calculate derivative factor - c = 1 - for _, i in zip(range(sum(delta)), range(sum(alpha), 0, -1)): - c *= i - alpha_ = numpy.array(alpha) - numpy.array(delta) - return c * bernstein_b(points, alpha_) - - -def bernstein_Db(points, alpha, order): +def bernstein_Dx(points, ks, order, R2B): points = numpy.asarray(points) - alpha = tuple(alpha) + ks = tuple(ks) N, d_1 = points.shape - assert d_1 == len(alpha) - Dshape = (d_1,) * order - - result = numpy.empty(Dshape + (N,)) - for indices in numpy.ndindex(Dshape): - delta = [0] * d_1 - for i in indices: - delta[i] += 1 - result[indices + (slice(None),)] = bernstein_db(points, alpha, delta) - return result - - -def bernstein_dx(points, alpha, delta, R2B): - points = numpy.asarray(points) - alpha = tuple(alpha) - delta = tuple(delta) - - N, d_1 = points.shape - assert d_1 == len(alpha) == len(delta) + 1 - - result = bernstein_Db(points, alpha, sum(delta)) - for d, c in enumerate(delta): - for _ in range(c): - result = R2B[:, d].dot(result) + 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) + Db_shape = (d_1,) * order + Db_tensor = numpy.empty(Db_shape + (N,)) + 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 From f8cafd0107ad938547e5a00bfa98bd3d9695df30 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Tue, 23 Oct 2018 15:03:02 +0200 Subject: [PATCH 09/46] reorder basis functions according to 'mis' --- FIAT/bernstein.py | 20 +++++++------------- test/unit/test_bernstein.py | 18 +++++++++--------- 2 files changed, 16 insertions(+), 22 deletions(-) diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index 39c7dbcc7..5b926a02c 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -17,9 +17,7 @@ # You should have received a copy of the GNU General Public License # along with FIAT. If not, see . -import itertools import math - import numpy from FIAT.finite_element import FiniteElement @@ -44,14 +42,12 @@ def __init__(self, ref_el, degree): # Generate triangular barycentric indices dim = ref_el.get_spatial_dimension() - alphas = [(degree - sum(beta),) + tuple(reversed(beta)) - for beta in itertools.product(range(degree + 1), repeat=dim) - if sum(beta) <= degree] + kss = mis(dim + 1, degree) # Fill data structures nodes = [] - for i, alpha in enumerate(alphas): - vertices, = numpy.nonzero(alpha) + 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) @@ -108,16 +104,14 @@ def tabulate(self, order, points, entity=None): # Generate triangular barycentric indices deg = self.degree() dim = ref_el.get_spatial_dimension() - etas = [(deg - sum(beta),) + tuple(reversed(beta)) - for beta in itertools.product(range(deg + 1), repeat=dim) - if sum(beta) <= deg] + kss = mis(dim + 1, deg) result = {} for D in range(order + 1): for alpha in mis(dim, D): - result[alpha] = numpy.zeros((len(etas), len(cell_points))) - for i, eta in enumerate(etas): - for alpha, vec in bernstein_Dx(B, eta, D, R2B).items(): + result[alpha] = numpy.zeros((len(kss), len(cell_points))) + for i, ks in enumerate(kss): + for alpha, vec in bernstein_Dx(B, ks, D, R2B).items(): result[alpha][i, :] = vec return result diff --git a/test/unit/test_bernstein.py b/test/unit/test_bernstein.py index 74de5d03c..301c427da 100644 --- a/test/unit/test_bernstein.py +++ b/test/unit/test_bernstein.py @@ -28,12 +28,12 @@ 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., 0., 0., 0., 0., 0.], - [0., 0., 0., 0., 0., 0.], [0.0831321, -2.12896637, 2.64569763, -7.25409741, 1.17096531, -6.51673126], - [-7.90833147, -7.90833147, -2.78320042, -2.78320042, -1.30846811, -1.30846811], [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], ]) @@ -41,12 +41,12 @@ 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], - [-3.95416573, -3.95416573, -1.39160021, -1.39160021, -0.65423405, -0.65423405], - [0., 0., 0., 0., 0., 0.], [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], - [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.], + [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.], ]) @@ -54,12 +54,12 @@ 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], - [-7.25409741, -6.51673126, -2.12896637, 1.17096531, 0.0831321, 2.64569763], - [3.95416573, 3.95416573, 1.39160021, 1.39160021, 0.65423405, 0.65423405], [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], - [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.], [0., 0., 0., 0., 0., 0.], ]) From d264e22f9449c19ab4ee77124536b8a8b51592cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Tue, 23 Oct 2018 15:43:48 +0200 Subject: [PATCH 10/46] add a docstring --- FIAT/bernstein.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index 5b926a02c..47331b1a0 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -152,6 +152,18 @@ def bernstein_db(points, ks, alpha=None): 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) From 90b49d276a9cb879209a0e2e6d4661ee1cedce19 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Tue, 23 Oct 2018 16:05:30 +0200 Subject: [PATCH 11/46] restructure to allow point evaluation of Bernstein elements --- FIAT/bernstein.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index 47331b1a0..8d0ef06c2 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -101,18 +101,22 @@ def tabulate(self, order, points, entity=None): B = numpy.hstack([cell_points, numpy.ones((len(cell_points), 1))]).dot(R2B.T) - # Generate triangular barycentric indices + # Evaluate everything deg = self.degree() dim = ref_el.get_spatial_dimension() - kss = mis(dim + 1, deg) - - result = {} - for D in range(order + 1): - for alpha in mis(dim, D): - result[alpha] = numpy.zeros((len(kss), len(cell_points))) - for i, ks in enumerate(kss): - for alpha, vec in bernstein_Dx(B, ks, D, R2B).items(): - result[alpha][i, :] = vec + 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 @@ -175,8 +179,9 @@ def bernstein_Dx(points, ks, order, R2B): 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,)) + Db_tensor = numpy.empty(Db_shape + (N,), dtype=dtype) for ds in numpy.ndindex(Db_shape): alpha = [0] * d_1 for d in ds: From c4b7b1b3a2089d72b95d9fd97b65860e239338da Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mikl=C3=B3s=20Homolya?= Date: Wed, 24 Oct 2018 08:23:58 +0200 Subject: [PATCH 12/46] add missing word 'Lesser' in license name --- FIAT/bernstein.py | 4 ++-- test/unit/test_bernstein.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/FIAT/bernstein.py b/FIAT/bernstein.py index 8d0ef06c2..251f94204 100644 --- a/FIAT/bernstein.py +++ b/FIAT/bernstein.py @@ -14,8 +14,8 @@ # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public # License for more details. # -# You should have received a copy of the GNU General Public License -# along with FIAT. If not, see . +# You should have received a copy of the GNU Lesser General Public +# License along with FIAT. If not, see . import math import numpy diff --git a/test/unit/test_bernstein.py b/test/unit/test_bernstein.py index 301c427da..2b4714a8a 100644 --- a/test/unit/test_bernstein.py +++ b/test/unit/test_bernstein.py @@ -14,8 +14,8 @@ # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public # License for more details. # -# You should have received a copy of the GNU General Public License -# along with FIAT. If not, see . +# You should have received a copy of the GNU Lesser General Public +# License along with FIAT. If not, see . import numpy import pytest From 8704e0b6b0554ff2daf93e8ba94cb3650281d1a0 Mon Sep 17 00:00:00 2001 From: Chris Richardson Date: Fri, 11 Jan 2019 19:13:54 +0000 Subject: [PATCH 13/46] Merged in chris/update-pytest (pull request #53) updates for latest pytest * Updates for pytest version * Update pipelines Approved-by: Garth Wells --- bitbucket-pipelines.yml | 4 +- test/regression/test_regression.py | 2 +- test/unit/test_fiat.py | 64 ++++++++++++++--------------- test/unit/test_quadrature.py | 40 +++++++++++++----- test/unit/test_reference_element.py | 18 ++++---- 5 files changed, 73 insertions(+), 55 deletions(-) diff --git a/bitbucket-pipelines.yml b/bitbucket-pipelines.yml index a32b26f1c..e3488e10e 100644 --- a/bitbucket-pipelines.yml +++ b/bitbucket-pipelines.yml @@ -6,9 +6,9 @@ pipelines: script: - apt-get update - apt-get install -y - git python3-minimal python3-pip python3-pytest python3-setuptools python3-wheel + git python3-minimal python3-pip python3-setuptools python3-wheel --no-install-recommends - - pip3 install flake8 + - pip3 install flake8 pytest --upgrade - pip3 install . - python3 -m flake8 - DATA_REPO_GIT="" python3 -m pytest -v test/ diff --git a/test/regression/test_regression.py b/test/regression/test_regression.py index f2f510f41..4d8afc572 100644 --- a/test/regression/test_regression.py +++ b/test/regression/test_regression.py @@ -277,7 +277,7 @@ def _perform_test(family, dim, degree, reference_table): reference = json.load(open(filename, "r"), object_hook=json_numpy_obj_hook) for test_case in test_cases: family, dim, degree = test_case - yield _perform_test, family, dim, degree, reference[str(test_case)] + _perform_test(family, dim, degree, reference[str(test_case)]) # Update references if missing except (IOError, KeyError) as e: diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index eec0bd817..e949d0468 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -203,42 +203,42 @@ def __init__(self, a, b): # MixedElement made of nodal elements should be nodal, but its API # is currently just broken. - xfail_impl("MixedElement([" - " DiscontinuousLagrange(T, 1)," - " RaviartThomas(T, 2)" - "])"), + pytest.param("MixedElement([" + " DiscontinuousLagrange(T, 1)," + " RaviartThomas(T, 2)" + "])", marks=xfail_impl), # Following element do not bother implementing get_nodal_basis # so the test would need to be rewritten using tabulate - xfail_impl("TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2))"), - xfail_impl("Hdiv(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"), - xfail_impl("Hcurl(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"), - xfail_impl("HDivTrace(T, 1)"), - xfail_impl("EnrichedElement(" - "Hdiv(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " - "Hdiv(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" - ")"), - xfail_impl("EnrichedElement(" - "Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " - "Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" - ")"), + pytest.param("TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2))", marks=xfail_impl), + pytest.param("Hdiv(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))", marks=xfail_impl), + pytest.param("Hcurl(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))", marks=xfail_impl), + pytest.param("HDivTrace(T, 1)", marks=xfail_impl), + pytest.param("EnrichedElement(" + "Hdiv(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " + "Hdiv(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" + ")", marks=xfail_impl), + pytest.param("EnrichedElement(" + "Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " + "Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" + ")", marks=xfail_impl), # Following elements are checked using tabulate - xfail_impl("HDivTrace(T, 0)"), - xfail_impl("HDivTrace(T, 1)"), - xfail_impl("HDivTrace(T, 2)"), - xfail_impl("HDivTrace(T, 3)"), - xfail_impl("HDivTrace(S, 0)"), - xfail_impl("HDivTrace(S, 1)"), - xfail_impl("HDivTrace(S, 2)"), - xfail_impl("HDivTrace(S, 3)"), - xfail_impl("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))"), - xfail_impl("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))"), - xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))"), - xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)), Lagrange(I, 2))"), - xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)))"), - xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)))"), - xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))"), - xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))"), + pytest.param("HDivTrace(T, 0)", marks=xfail_impl), + pytest.param("HDivTrace(T, 1)", marks=xfail_impl), + pytest.param("HDivTrace(T, 2)", marks=xfail_impl), + pytest.param("HDivTrace(T, 3)", marks=xfail_impl), + pytest.param("HDivTrace(S, 0)", marks=xfail_impl), + pytest.param("HDivTrace(S, 1)", marks=xfail_impl), + pytest.param("HDivTrace(S, 2)", marks=xfail_impl), + pytest.param("HDivTrace(S, 3)", marks=xfail_impl), + pytest.param("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))", marks=xfail_impl), + pytest.param("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))", marks=xfail_impl), + pytest.param("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))", marks=xfail_impl), + pytest.param("TensorProductElement(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)), Lagrange(I, 2))", marks=xfail_impl), + pytest.param("FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)))", marks=xfail_impl), + pytest.param("FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)))", marks=xfail_impl), + pytest.param("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))", marks=xfail_impl), + pytest.param("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))", marks=xfail_impl), ] diff --git a/test/unit/test_quadrature.py b/test/unit/test_quadrature.py index 68c5c7e1c..fdcf88e66 100644 --- a/test/unit/test_quadrature.py +++ b/test/unit/test_quadrature.py @@ -134,21 +134,39 @@ def test_create_quadrature_extr_quadrilateral(extr_quadrilateral, basedeg, extrd (2**(basedeg + 2) - 2) / ((basedeg + 1)*(basedeg + 2)) * 1/(extrdeg + 1)) -@pytest.mark.parametrize("cell", [interval(), - triangle(), - tetrahedron(), - quadrilateral()]) -def test_invalid_quadrature_degree(cell, scheme): +def test_invalid_quadrature_degree_interval(interval, scheme): with pytest.raises(ValueError): - FIAT.create_quadrature(cell, -1, scheme) + FIAT.create_quadrature(interval, -1, scheme) -@pytest.mark.parametrize("cell", [extr_interval(), - extr_triangle(), - extr_quadrilateral()]) -def test_invalid_quadrature_degree_tensor_prod(cell): +def test_invalid_quadrature_degree_tri(triangle, scheme): with pytest.raises(ValueError): - FIAT.create_quadrature(cell, (-1, -1)) + FIAT.create_quadrature(triangle, -1, scheme) + + +def test_invalid_quadrature_degree_quad(quadrilateral, scheme): + with pytest.raises(ValueError): + FIAT.create_quadrature(quadrilateral, -1, scheme) + + +def test_invalid_quadrature_degree_tet(tetrahedron, scheme): + with pytest.raises(ValueError): + FIAT.create_quadrature(tetrahedron, -1, scheme) + + +def test_invalid_quadrature_degree_tensor_prod_interval(extr_interval): + with pytest.raises(ValueError): + FIAT.create_quadrature(extr_interval, (-1, -1)) + + +def test_invalid_quadrature_degree_tensor_prod_tri(extr_triangle): + with pytest.raises(ValueError): + FIAT.create_quadrature(extr_triangle, (-1, -1)) + + +def test_invalid_quadrature_degree_tensor_prod_quad(extr_quadrilateral): + with pytest.raises(ValueError): + FIAT.create_quadrature(extr_quadrilateral, (-1, -1)) def test_tensor_product_composition(interval, triangle, extr_triangle, scheme): diff --git a/test/unit/test_reference_element.py b/test/unit/test_reference_element.py index edc1fea65..a509b6c64 100644 --- a/test/unit/test_reference_element.py +++ b/test/unit/test_reference_element.py @@ -39,8 +39,8 @@ @pytest.mark.parametrize(('cell', 'connectivity'), [(tetrahedron, ufc_tetrahedron_21connectivity), (hexahedron, ufc_hexahedron_21connectivity), - pytest.mark.xfail((triangle_x_interval, [])), - pytest.mark.xfail((quadrilateral_x_interval, []))]) + pytest.param(triangle_x_interval, [], marks=pytest.mark.xfail), + pytest.param(quadrilateral_x_interval, [], marks=pytest.mark.xfail)]) def test_ufc_connectivity_21(cell, connectivity): """Check face-edge connectivity builds what UFC expects. This is only non-trivial case ; the rest is x-0 and D-x, @@ -51,9 +51,9 @@ def test_ufc_connectivity_21(cell, connectivity): @pytest.mark.parametrize('cell', [point, interval, triangle, tetrahedron, quadrilateral, hexahedron, - pytest.mark.xfail(interval_x_interval), - pytest.mark.xfail(triangle_x_interval), - pytest.mark.xfail(quadrilateral_x_interval)]) + pytest.param(interval_x_interval, marks=pytest.mark.xfail), + pytest.param(triangle_x_interval, marks=pytest.mark.xfail), + pytest.param(quadrilateral_x_interval, marks=pytest.mark.xfail)]) def test_ufc_connectivity_x0(cell): """Check x-0 connectivity is just what get_topology gives""" for dim0 in range(cell.get_spatial_dimension()+1): @@ -66,9 +66,9 @@ def test_ufc_connectivity_x0(cell): @pytest.mark.parametrize('cell', [point, interval, triangle, tetrahedron, quadrilateral, hexahedron, - pytest.mark.xfail(interval_x_interval), - pytest.mark.xfail(triangle_x_interval), - pytest.mark.xfail(quadrilateral_x_interval)]) + pytest.param(interval_x_interval, marks=pytest.mark.xfail), + pytest.param(triangle_x_interval, marks=pytest.mark.xfail), + pytest.param(quadrilateral_x_interval, marks=pytest.mark.xfail)]) def test_ufc_connectivity_Dx(cell): """Check D-x connectivity is just [(0,1,2,...)]""" D = cell.get_spatial_dimension() @@ -79,7 +79,7 @@ def test_ufc_connectivity_Dx(cell): @pytest.mark.parametrize(('cell', 'volume'), - [pytest.mark.xfail((point, 1)), + [pytest.param(point, 1, marks=pytest.mark.xfail), (interval, 1), (triangle, 1/2), (quadrilateral, 1), From 102128721a18422adeaa8e242676fb7dd656e464 Mon Sep 17 00:00:00 2001 From: Cyrus Cheng Date: Mon, 28 Jan 2019 13:44:00 +0000 Subject: [PATCH 14/46] Add new discontinuous lagrange cube module. --- FIAT/discontinuous_lagrange_cube.py | 55 +++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) create mode 100644 FIAT/discontinuous_lagrange_cube.py diff --git a/FIAT/discontinuous_lagrange_cube.py b/FIAT/discontinuous_lagrange_cube.py new file mode 100644 index 000000000..a7286e083 --- /dev/null +++ b/FIAT/discontinuous_lagrange_cube.py @@ -0,0 +1,55 @@ +# 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 . + +from FIAT import finite_element, polynomial_set, dual_set, functional, P0, + reference_element +from FIAT.discontinuous_lagrange import DiscontinuousLagrangeDualSet + + +class DiscontinuousLagrangeCubeDualSet(DiscontinuousLagrangeDualSet): + """The dual basis for Lagrange 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): + super(DiscontinuousLagrangeCubeDualSet, self).__init__(ref_el, degree) + + +class HigherOrderDiscontinuousLagrangeCube(finite_element.CiarletElement): + """The discontinuous Lagrange finite element. It is what it is.""" + + def __init__(self, ref_el, degree): + hypercube_simplex_map = {Point: Point, + DefaultLine: DefaultLine, + UFCInterval: UFCInterval, + UFCQuadrilateral: UFCTriangle, + UFCHexahedron: UFCTetrahedron} + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], degree) + dual = DiscontinuousLagrangeCubeDualSet(ref_el, degree) + formdegree = ref_el.get_spatial_dimension() # n-form + super(HigherOrderDiscontinuousLagrangeCube, self).__init__(poly_set, dual, degree, formdegree) + + +def DiscontinuousLagrangeCube(ref_el, degree): + if degree == 0: + hypercube_simplex_map = {Point: Point, + DefaultLine: DefaultLine, + UFCInterval: UFCInterval, + UFCQuadrilateral: UFCTriangle, + UFCHexahedron: UFCTetrahedron} + return P0.P0(hypercube_simplex_map[ref_el]) + else: + return HigherOrderDiscontinuousLagrangeCube(ref_el, degree) From 55fea3ded6ec52bf349d519e4eebae38e7a972f3 Mon Sep 17 00:00:00 2001 From: Cyrus Cheng Date: Mon, 28 Jan 2019 14:37:57 +0000 Subject: [PATCH 15/46] Add ref_el as optional argument to CiarletElement. --- FIAT/finite_element.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index 7ca937534..aba117298 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, ref_el=None, formdegree=None, mapping="affine"): + ref_el = ref_el or poly_set.get_reference_element() super(CiarletElement, self).__init__(ref_el, dual, order, formdegree, mapping) # build generalized Vandermonde matrix From ba9f4a567289814e975373f26d670ec421e0e9f0 Mon Sep 17 00:00:00 2001 From: Cyrus Cheng Date: Mon, 28 Jan 2019 14:41:35 +0000 Subject: [PATCH 16/46] Renamed discontinuous serendipity element. --- FIAT/discontinuous_lagrange_cube.py | 55 ---------------------- FIAT/discontinuous_serendipity.py | 72 +++++++++++++++++++++++++++++ 2 files changed, 72 insertions(+), 55 deletions(-) delete mode 100644 FIAT/discontinuous_lagrange_cube.py create mode 100644 FIAT/discontinuous_serendipity.py diff --git a/FIAT/discontinuous_lagrange_cube.py b/FIAT/discontinuous_lagrange_cube.py deleted file mode 100644 index a7286e083..000000000 --- a/FIAT/discontinuous_lagrange_cube.py +++ /dev/null @@ -1,55 +0,0 @@ -# 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 . - -from FIAT import finite_element, polynomial_set, dual_set, functional, P0, - reference_element -from FIAT.discontinuous_lagrange import DiscontinuousLagrangeDualSet - - -class DiscontinuousLagrangeCubeDualSet(DiscontinuousLagrangeDualSet): - """The dual basis for Lagrange 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): - super(DiscontinuousLagrangeCubeDualSet, self).__init__(ref_el, degree) - - -class HigherOrderDiscontinuousLagrangeCube(finite_element.CiarletElement): - """The discontinuous Lagrange finite element. It is what it is.""" - - def __init__(self, ref_el, degree): - hypercube_simplex_map = {Point: Point, - DefaultLine: DefaultLine, - UFCInterval: UFCInterval, - UFCQuadrilateral: UFCTriangle, - UFCHexahedron: UFCTetrahedron} - poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], degree) - dual = DiscontinuousLagrangeCubeDualSet(ref_el, degree) - formdegree = ref_el.get_spatial_dimension() # n-form - super(HigherOrderDiscontinuousLagrangeCube, self).__init__(poly_set, dual, degree, formdegree) - - -def DiscontinuousLagrangeCube(ref_el, degree): - if degree == 0: - hypercube_simplex_map = {Point: Point, - DefaultLine: DefaultLine, - UFCInterval: UFCInterval, - UFCQuadrilateral: UFCTriangle, - UFCHexahedron: UFCTetrahedron} - return P0.P0(hypercube_simplex_map[ref_el]) - else: - return HigherOrderDiscontinuousLagrangeCube(ref_el, degree) diff --git a/FIAT/discontinuous_serendipity.py b/FIAT/discontinuous_serendipity.py new file mode 100644 index 000000000..1370a99be --- /dev/null +++ b/FIAT/discontinuous_serendipity.py @@ -0,0 +1,72 @@ +# 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 . + +from FIAT import finite_element, polynomial_set, dual_set, functional, P0, + reference_element +from FIAT.discontinuous_lagrange import DiscontinuousLagrangeDualSet +from FIAT.P0 import P0Dual + +hypercube_simplex_map = {Point: Point, + DefaultLine: DefaultLine, + UFCInterval: UFCInterval, + UFCQuadrilateral: UFCTriangle, + UFCHexahedron: UFCTetrahedron} + +class DPC0Dual(P0Dual): + def __init__(self, ref_el): + super(DPC0Dual, self).__init__(ref_el) + +class DPC0(finite_element.CiarletElement): + def __init__(self, ref_el): + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) + dual = DPC0Dual(ref_el) + degree = 0 + formdegree = ref_el.get_spatial_dimension() # n-form + super(DPC0, self).__init__(poly_set=poly_set, + dual=dual, + degree=degree, + ref_el=ref_el, + formdegree=formdegree) + + +class DiscontinuousSerendipityDualSet(DiscontinuousLagrangeDualSet): + """The dual basis for Lagrange 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): + super(DiscontinuousSerendipityDualSet, self).__init__(ref_el, degree) + + +class HigherOrderDiscontinuousSerendipity(finite_element.CiarletElement): + """The discontinuous Lagrange 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 = DiscontinuousSerendipityDualSet(hypercube_simplex_map[ref_el], degree) + formdegree = ref_el.get_spatial_dimension() # n-form + super(HigherOrderDiscontinuousSerendipity, self).__init__(poly_set=poly_set, + dual=dual, + degree=degree, + ref_el=ref_el, + formdegree=formdegree) + + +def DiscontinuousSerendipity(ref_el, degree): + if degree == 0: + return DPC0(ref_el) + else: + return HigherOrderDiscontinuousSerendipity(ref_el, degree) From 4f73c356c25811cbf63826b6790d80de5da02206 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Mon, 28 Jan 2019 19:17:40 +0000 Subject: [PATCH 17/46] Made changes to hypercube_simplex_map. --- FIAT/discontinuous_serendipity.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/FIAT/discontinuous_serendipity.py b/FIAT/discontinuous_serendipity.py index 1370a99be..0ec8377fc 100644 --- a/FIAT/discontinuous_serendipity.py +++ b/FIAT/discontinuous_serendipity.py @@ -13,16 +13,16 @@ # You should have received a copy of the GNU Lesser General Public License # along with FIAT. If not, see . -from FIAT import finite_element, polynomial_set, dual_set, functional, P0, - reference_element +from FIAT import finite_element, polynomial_set, dual_set, functional, P0, reference_element from FIAT.discontinuous_lagrange import DiscontinuousLagrangeDualSet +from FIAT.reference_element import Point, DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron, UFCTriangle, UFCTetrahedron from FIAT.P0 import P0Dual -hypercube_simplex_map = {Point: Point, - DefaultLine: DefaultLine, - UFCInterval: UFCInterval, - UFCQuadrilateral: UFCTriangle, - UFCHexahedron: UFCTetrahedron} +hypercube_simplex_map = {'Point': Point(), + 'DefaultLine': DefaultLine(), + 'UFCInterval': UFCInterval(), + 'UFCQuadrilateral': UFCTriangle(), + 'UFCHexahedron': UFCTetrahedron()} class DPC0Dual(P0Dual): def __init__(self, ref_el): @@ -30,19 +30,19 @@ def __init__(self, ref_el): class DPC0(finite_element.CiarletElement): def __init__(self, ref_el): - poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[type(ref_el).__name__], 0) dual = DPC0Dual(ref_el) degree = 0 formdegree = ref_el.get_spatial_dimension() # n-form super(DPC0, self).__init__(poly_set=poly_set, dual=dual, - degree=degree, + order=degree, ref_el=ref_el, formdegree=formdegree) class DiscontinuousSerendipityDualSet(DiscontinuousLagrangeDualSet): - """The dual basis for Lagrange elements. This class works for + """The dual basis for Serendipity 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""" @@ -52,15 +52,15 @@ def __init__(self, ref_el, degree): class HigherOrderDiscontinuousSerendipity(finite_element.CiarletElement): - """The discontinuous Lagrange finite element. It is what it is.""" + """The discontinuous Serendipity 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 = DiscontinuousSerendipityDualSet(hypercube_simplex_map[ref_el], degree) + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[type(ref_el).__name__], degree) + dual = DiscontinuousSerendipityDualSet(hypercube_simplex_map[type(ref_el).__name__], degree) formdegree = ref_el.get_spatial_dimension() # n-form super(HigherOrderDiscontinuousSerendipity, self).__init__(poly_set=poly_set, dual=dual, - degree=degree, + order=degree, ref_el=ref_el, formdegree=formdegree) From f5a765302e732e10eda56cadaba7b4eb0c30b6ab Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Thu, 3 Jan 2019 13:59:50 -0600 Subject: [PATCH 18/46] general/abstract to_riesz, lift to dual --- FIAT/dual_set.py | 77 +++++++++++++++++++++++++++++- FIAT/functional.py | 116 +++++++++++++++++++++++++++++++++++++-------- 2 files changed, 172 insertions(+), 21 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index b9e2e74e2..f24d7473d 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -15,7 +15,9 @@ # You should have received a copy of the GNU Lesser General Public License # along with FIAT. If not, see . -import numpy +import numpy, collections + +from FIAT import polynomial_set class DualSet(object): @@ -64,3 +66,76 @@ def to_riesz(self, poly_set): self.mat[i][:] = self.nodes[i].to_riesz(poly_set) return self.mat + + def to_riesz_foo(self, poly_set): + + tshape = self.nodes[0].target_shape + num_nodes = len(self.nodes) + es = poly_set.get_expansion_set() + ed = poly_set.get_embedded_degree() + num_exp = es.get_num_members(poly_set.get_embedded_degree()) + nbf = poly_set.get_num_members() + + riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp]) + + self.mat = numpy.zeros(riesz_shape, "d") + + # let's amalgamate the pt_dict and deriv_dicts of all the + # functionals so we can tabulate the basis functions twice only + # (once on pts and once on derivatives) + + # Need: dictionary mapping pts to which functionals they come from + pts_to_ells = collections.OrderedDict() + dpts_to_ells = collections.OrderedDict() + + for i, ell in enumerate(self.nodes): + for pt in ell.pt_dict: + if pt in pts_to_ells: + pts_to_ells[pt].append(i) + else: + pts_to_ells[pt] = [i] + + for i, ell in enumerate(self.nodes): + for pt in ell.deriv_dict: + if pt in dpts_to_ells: + dpts_to_ells[pt].append(i) + else: + dpts_to_ells[pt] = [i] + + # Now tabulate + pts = list(pts_to_ells.keys()) + bfs = es.tabulate(ed, pts) + + + for j, pt in enumerate(pts): + which_ells = pts_to_ells[pt] + + for k in which_ells: + pt_dict = self.nodes[k].pt_dict + wc_list = pt_dict[pt] + + for i in range(bfs.shape[0]): + for (w, c) in wc_list: + + self.mat[k][c][i] += w*bfs[i, j] + + + mdo = max([ell.max_deriv_order for ell in self.nodes]) + if mdo > 0: + dpts = list(dpts_to_ells.keys()) + es_foo = polynomial_set.ONPolynomialSet(self.ref_el, ed) + dbfs = es_foo.tabulate(dpts, mdo) + + for j, pt in enumerate(dpts): + which_ells = dpts_to_ells[pt] + + for k in which_ells: + dpt_dict = self.nodes[k].deriv_dict + wac_list = dpt_dict[pt] + + for i in range(nbf): + for (w, alpha, c) in wac_list: + self.mat[k][c][i] += w*dbfs[alpha][i, j] + + #assert numpy.allclose(self.mat, self.to_riesz_old(poly_set)) + return self.mat diff --git a/FIAT/functional.py b/FIAT/functional.py index 8d66327a2..441457c1e 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -26,6 +26,8 @@ import numpy import sympy +from FIAT import polynomial_set + def index_iterator(shp): """Constructs a generator iterating over all indices in @@ -42,6 +44,7 @@ def index_iterator(shp): for foo in index_iterator(shp_foo): yield [i] + foo + # also put in a "jet_dict" that maps # pt --> {wt, multiindex, comp} # the multiindex is an iterable of nonnegative @@ -103,17 +106,20 @@ def to_riesz(self, poly_set): ed = poly_set.get_embedded_degree() pt_dict = self.get_point_dict() + nbf = poly_set.get_num_members() + pts = list(pt_dict.keys()) # bfs is matrix that is pdim rows by num_pts cols # where pdim is the polynomial dimension bfs = es.tabulate(ed, pts) + npts = len(pts) result = numpy.zeros(poly_set.coeffs.shape[1:], "d") # loop over points - for j in range(len(pts)): + for j in range(npts): pt_cur = pts[j] wc_list = pt_dict[pt_cur] @@ -123,7 +129,25 @@ def to_riesz(self, poly_set): result[c][i] += w * bfs[i, j] if self.deriv_dict: - raise NotImplementedError("Generic to_riesz implementation does not support derivatives") + dpt_dict = self.deriv_dict + + # this makes things quicker since it uses dmats after + # instantiation + es_foo = polynomial_set.ONPolynomialSet(self.ref_el, ed) + dpts = list(dpt_dict.keys()) + + # figure out the maximum order of derivative being taken + dbfs = es_foo.tabulate(dpts, self.max_deriv_order) + + ndpts = len(dpts) + for j in range(ndpts): + dpt_cur = dpts[j] + wac_list = dpt_dict[dpt_cur] + for i in range(nbf): + for (w, alpha, c) in wac_list: + result[c][i] += w * dbfs[alpha][i, j] + + #raise NotImplementedError("Generic to_riesz implementation does not support derivatives") return result @@ -172,8 +196,8 @@ class PointDerivative(Functional): functions at a particular point x.""" def __init__(self, ref_el, x, alpha): - dpt_dict = {x: [(1.0, alpha, tuple())]} - self.alpha = alpha + dpt_dict = {x: [(1.0, tuple(alpha), tuple())]} + self.alpha = tuple(alpha) self.order = sum(self.alpha) Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointDeriv") @@ -191,24 +215,24 @@ def __call__(self, fn): return sympy.diff(fn(X), *dvars).evalf(subs=dict(zip(dX, x))) - def to_riesz(self, poly_set): - x = list(self.deriv_dict.keys())[0] + # def to_riesz(self, poly_set): + # x = list(self.deriv_dict.keys())[0] - X = sympy.DeferredVector('x') - dx = numpy.asarray([X[i] for i in range(len(x))]) + # X = sympy.DeferredVector('x') + # dx = numpy.asarray([X[i] for i in range(len(x))]) - es = poly_set.get_expansion_set() - ed = poly_set.get_embedded_degree() + # es = poly_set.get_expansion_set() + # ed = poly_set.get_embedded_degree() - bfs = es.tabulate(ed, [dx])[:, 0] + # bfs = es.tabulate(ed, [dx])[:, 0] - # Expand the multi-index as a series of variables to - # differentiate with respect to. - dvars = tuple(d for d, a in zip(dx, self.alpha) - for count in range(a)) + # # Expand the multi-index as a series of variables to + # # differentiate with respect to. + # dvars = tuple(d for d, a in zip(dx, self.alpha) + # for count in range(a)) - return numpy.asarray([sympy.lambdify(X, sympy.diff(b, *dvars))(x) - for b in bfs]) + # return numpy.asarray([sympy.lambdify(X, sympy.diff(b, *dvars))(x) + # for b in bfs]) class PointNormalDerivative(Functional): @@ -223,11 +247,55 @@ def __init__(self, ref_el, facet_no, pt): alpha = [0] * sd alpha[i] = 1 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") + + # def to_riesz(self, poly_set): + # x = list(self.deriv_dict.keys())[0] + + # X = sympy.DeferredVector('x') + # dx = numpy.asarray([X[i] for i in range(len(x))]) + + # es = poly_set.get_expansion_set() + # ed = poly_set.get_embedded_degree() + + # bfs = es.tabulate(ed, [dx])[:, 0] + + # # We need the gradient dotted with the normal. + # return numpy.asarray( + # [sympy.lambdify( + # X, sum([sympy.diff(b, dxi)*ni + # for dxi, ni in zip(dx, self.n)]))(x) + # for b in bfs]) + + +class PointNormalSecondDerivative(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() + tau = numpy.zeros((sd*(sd+1)//2,)) + + alphas = [] + cur = 0 + for i in range(sd): + for j in range(i, sd): + alpha = [0] * sd + alpha[i] += 1 + alpha[j] += 1 + alphas.append(tuple(alpha)) + tau[cur] = n[i]*n[j] + cur += 1 + + self.tau = tau + 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") - def to_riesz(self, poly_set): + def to_riesz_blah(self, poly_set): x = list(self.deriv_dict.keys())[0] X = sympy.DeferredVector('x') @@ -239,10 +307,18 @@ def to_riesz(self, poly_set): bfs = es.tabulate(ed, [dx])[:, 0] # We need the gradient dotted with the normal. + def diffargtosplat(alpha): + assert sum(alpha) > 0 + args = [] + for dxi, a in zip(dx, alpha): + if a > 0: + args.extend([dxi, a]) + return tuple(args) + return numpy.asarray( [sympy.lambdify( - X, sum([sympy.diff(b, dxi)*ni - for dxi, ni in zip(dx, self.n)]))(x) + X, sum([sympy.diff(b, *diffargtosplat(alpha))*taui + for alpha, taui in zip(self.alphas, self.tau)]))(x) for b in bfs]) From d9ddd6039caaefee121794ca60ca95b0eeffb12e Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Mon, 14 Jan 2019 13:05:00 -0600 Subject: [PATCH 19/46] Switch to new to_riesz method --- FIAT/dual_set.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index f24d7473d..d45c9501b 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -51,23 +51,23 @@ def get_entity_ids(self): def get_reference_element(self): return self.ref_el - def to_riesz(self, poly_set): + # def to_riesz(self, poly_set): - tshape = self.nodes[0].target_shape - num_nodes = len(self.nodes) - es = poly_set.get_expansion_set() - num_exp = es.get_num_members(poly_set.get_embedded_degree()) + # tshape = self.nodes[0].target_shape + # num_nodes = len(self.nodes) + # es = poly_set.get_expansion_set() + # num_exp = es.get_num_members(poly_set.get_embedded_degree()) - riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp]) + # riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp]) - self.mat = numpy.zeros(riesz_shape, "d") + # self.mat = numpy.zeros(riesz_shape, "d") - for i in range(len(self.nodes)): - self.mat[i][:] = self.nodes[i].to_riesz(poly_set) + # for i in range(len(self.nodes)): + # self.mat[i][:] = self.nodes[i].to_riesz(poly_set) - return self.mat + # return self.mat - def to_riesz_foo(self, poly_set): + def to_riesz(self, poly_set): tshape = self.nodes[0].target_shape num_nodes = len(self.nodes) From a80ebbdabc9f566c994d76b73292a59ba15fa631 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Mon, 14 Jan 2019 13:11:05 -0600 Subject: [PATCH 20/46] pep 8, remove some old code --- FIAT/dual_set.py | 22 ++-------------------- FIAT/functional.py | 4 +--- 2 files changed, 3 insertions(+), 23 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index d45c9501b..3977518f2 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -15,7 +15,8 @@ # You should have received a copy of the GNU Lesser General Public License # along with FIAT. If not, see . -import numpy, collections +import numpy +import collections from FIAT import polynomial_set @@ -51,22 +52,6 @@ def get_entity_ids(self): def get_reference_element(self): return self.ref_el - # def to_riesz(self, poly_set): - - # tshape = self.nodes[0].target_shape - # num_nodes = len(self.nodes) - # es = poly_set.get_expansion_set() - # num_exp = es.get_num_members(poly_set.get_embedded_degree()) - - # riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp]) - - # self.mat = numpy.zeros(riesz_shape, "d") - - # for i in range(len(self.nodes)): - # self.mat[i][:] = self.nodes[i].to_riesz(poly_set) - - # return self.mat - def to_riesz(self, poly_set): tshape = self.nodes[0].target_shape @@ -106,7 +91,6 @@ def to_riesz(self, poly_set): pts = list(pts_to_ells.keys()) bfs = es.tabulate(ed, pts) - for j, pt in enumerate(pts): which_ells = pts_to_ells[pt] @@ -119,7 +103,6 @@ def to_riesz(self, poly_set): self.mat[k][c][i] += w*bfs[i, j] - mdo = max([ell.max_deriv_order for ell in self.nodes]) if mdo > 0: dpts = list(dpts_to_ells.keys()) @@ -137,5 +120,4 @@ def to_riesz(self, poly_set): for (w, alpha, c) in wac_list: self.mat[k][c][i] += w*dbfs[alpha][i, j] - #assert numpy.allclose(self.mat, self.to_riesz_old(poly_set)) return self.mat diff --git a/FIAT/functional.py b/FIAT/functional.py index 441457c1e..43fa9750d 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -133,7 +133,7 @@ def to_riesz(self, poly_set): # this makes things quicker since it uses dmats after # instantiation - es_foo = polynomial_set.ONPolynomialSet(self.ref_el, ed) + es_foo = polynomial_set.ONPolynomialSet(self.ref_el, ed) dpts = list(dpt_dict.keys()) # figure out the maximum order of derivative being taken @@ -147,8 +147,6 @@ def to_riesz(self, poly_set): for (w, alpha, c) in wac_list: result[c][i] += w * dbfs[alpha][i, j] - #raise NotImplementedError("Generic to_riesz implementation does not support derivatives") - return result def tostr(self): From b2b5c9fe5b17d034feb8ba47cde6e6158c55ca04 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Tue, 15 Jan 2019 12:08:50 -0600 Subject: [PATCH 21/46] remove some code, comments --- FIAT/functional.py | 70 ---------------------------------------------- 1 file changed, 70 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index 43fa9750d..6eb527a5e 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -45,12 +45,6 @@ def index_iterator(shp): yield [i] + foo -# also put in a "jet_dict" that maps -# pt --> {wt, multiindex, comp} -# the multiindex is an iterable of nonnegative -# integers - - class Functional(object): """Class implementing an abstract functional. All functionals are discrete in the sense that @@ -97,7 +91,6 @@ def get_type_tag(self): normal component, which is probably handy for clients of FIAT""" return self.functional_type - # overload me in subclasses to make life easier!! def to_riesz(self, poly_set): """Constructs an array representation of the functional over the base of the given polynomial_set so that f(phi) for any @@ -213,25 +206,6 @@ def __call__(self, fn): return sympy.diff(fn(X), *dvars).evalf(subs=dict(zip(dX, x))) - # def to_riesz(self, poly_set): - # x = list(self.deriv_dict.keys())[0] - - # X = sympy.DeferredVector('x') - # dx = numpy.asarray([X[i] for i in range(len(x))]) - - # es = poly_set.get_expansion_set() - # ed = poly_set.get_embedded_degree() - - # bfs = es.tabulate(ed, [dx])[:, 0] - - # # Expand the multi-index as a series of variables to - # # differentiate with respect to. - # dvars = tuple(d for d, a in zip(dx, self.alpha) - # for count in range(a)) - - # return numpy.asarray([sympy.lambdify(X, sympy.diff(b, *dvars))(x) - # for b in bfs]) - class PointNormalDerivative(Functional): @@ -249,24 +223,6 @@ def __init__(self, ref_el, facet_no, pt): Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv") - # def to_riesz(self, poly_set): - # x = list(self.deriv_dict.keys())[0] - - # X = sympy.DeferredVector('x') - # dx = numpy.asarray([X[i] for i in range(len(x))]) - - # es = poly_set.get_expansion_set() - # ed = poly_set.get_embedded_degree() - - # bfs = es.tabulate(ed, [dx])[:, 0] - - # # We need the gradient dotted with the normal. - # return numpy.asarray( - # [sympy.lambdify( - # X, sum([sympy.diff(b, dxi)*ni - # for dxi, ni in zip(dx, self.n)]))(x) - # for b in bfs]) - class PointNormalSecondDerivative(Functional): @@ -293,32 +249,6 @@ def __init__(self, ref_el, facet_no, pt): Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv") - def to_riesz_blah(self, poly_set): - x = list(self.deriv_dict.keys())[0] - - X = sympy.DeferredVector('x') - dx = numpy.asarray([X[i] for i in range(len(x))]) - - es = poly_set.get_expansion_set() - ed = poly_set.get_embedded_degree() - - bfs = es.tabulate(ed, [dx])[:, 0] - - # We need the gradient dotted with the normal. - def diffargtosplat(alpha): - assert sum(alpha) > 0 - args = [] - for dxi, a in zip(dx, alpha): - if a > 0: - args.extend([dxi, a]) - return tuple(args) - - return numpy.asarray( - [sympy.lambdify( - X, sum([sympy.diff(b, *diffargtosplat(alpha))*taui - for alpha, taui in zip(self.alphas, self.tau)]))(x) - for b in bfs]) - class IntegralMoment(Functional): """An IntegralMoment is a functional""" From 7bd464b6ac7dd459f59e9f8d59496473b8885e5b Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Wed, 16 Jan 2019 11:22:40 -0600 Subject: [PATCH 22/46] Address dham's concerns on dual_set.py --- FIAT/dual_set.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index 3977518f2..fe77bc0a1 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -63,7 +63,7 @@ def to_riesz(self, poly_set): riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp]) - self.mat = numpy.zeros(riesz_shape, "d") + self.matrix = numpy.zeros(riesz_shape, "d") # let's amalgamate the pt_dict and deriv_dicts of all the # functionals so we can tabulate the basis functions twice only @@ -89,7 +89,7 @@ def to_riesz(self, poly_set): # Now tabulate pts = list(pts_to_ells.keys()) - bfs = es.tabulate(ed, pts) + expansion_values = es.tabulate(ed, pts) for j, pt in enumerate(pts): which_ells = pts_to_ells[pt] @@ -98,16 +98,19 @@ def to_riesz(self, poly_set): pt_dict = self.nodes[k].pt_dict wc_list = pt_dict[pt] - for i in range(bfs.shape[0]): + for i in range(expanion_values.shape[0]): for (w, c) in wc_list: - self.mat[k][c][i] += w*bfs[i, j] + self.matrix[k][c][i] += w*expanion_values[i, j] - mdo = max([ell.max_deriv_order for ell in self.nodes]) - if mdo > 0: + max_deriv_order = max([ell.max_deriv_order for ell in self.nodes]) + if max_deriv_order > 0: dpts = list(dpts_to_ells.keys()) - es_foo = polynomial_set.ONPolynomialSet(self.ref_el, ed) - dbfs = es_foo.tabulate(dpts, mdo) + # It's easiest/most efficient to get derivatives of the + # expansion set through the polynomial set interface. + # This is creating a short-lived set to do just this. + expansion = polynomial_set.ONPolynomialSet(self.ref_el, ed) + dexpanion_values = expansion.tabulate(dpts, max_deriv_order) for j, pt in enumerate(dpts): which_ells = dpts_to_ells[pt] @@ -118,6 +121,6 @@ def to_riesz(self, poly_set): for i in range(nbf): for (w, alpha, c) in wac_list: - self.mat[k][c][i] += w*dbfs[alpha][i, j] + self.matrix[k][c][i] += w*dexpanion_values[alpha][i, j] - return self.mat + return self.matrix From d94a88371844e2f9faf436c5f17f8fca40d8532a Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 18 Jan 2019 14:04:47 -0600 Subject: [PATCH 23/46] fix typo --- FIAT/dual_set.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index fe77bc0a1..e7e776161 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -98,10 +98,10 @@ def to_riesz(self, poly_set): pt_dict = self.nodes[k].pt_dict wc_list = pt_dict[pt] - for i in range(expanion_values.shape[0]): + for i in range(expansion_values.shape[0]): for (w, c) in wc_list: - self.matrix[k][c][i] += w*expanion_values[i, j] + self.matrix[k][c][i] += w*expansion_values[i, j] max_deriv_order = max([ell.max_deriv_order for ell in self.nodes]) if max_deriv_order > 0: @@ -110,7 +110,7 @@ def to_riesz(self, poly_set): # expansion set through the polynomial set interface. # This is creating a short-lived set to do just this. expansion = polynomial_set.ONPolynomialSet(self.ref_el, ed) - dexpanion_values = expansion.tabulate(dpts, max_deriv_order) + dexpansion_values = expansion.tabulate(dpts, max_deriv_order) for j, pt in enumerate(dpts): which_ells = dpts_to_ells[pt] @@ -121,6 +121,6 @@ def to_riesz(self, poly_set): for i in range(nbf): for (w, alpha, c) in wac_list: - self.matrix[k][c][i] += w*dexpanion_values[alpha][i, j] + self.matrix[k][c][i] += w*dexpansion_values[alpha][i, j] return self.matrix From 703b859c483c7aaae279c5fa7fe14da75bbc5184 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 18 Jan 2019 23:48:11 -0600 Subject: [PATCH 24/46] bug fix in derivatives for to_riesz --- FIAT/dual_set.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index e7e776161..276a112ee 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -60,7 +60,7 @@ def to_riesz(self, poly_set): ed = poly_set.get_embedded_degree() num_exp = es.get_num_members(poly_set.get_embedded_degree()) nbf = poly_set.get_num_members() - + riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp]) self.matrix = numpy.zeros(riesz_shape, "d") @@ -100,7 +100,6 @@ def to_riesz(self, poly_set): for i in range(expansion_values.shape[0]): for (w, c) in wc_list: - self.matrix[k][c][i] += w*expansion_values[i, j] max_deriv_order = max([ell.max_deriv_order for ell in self.nodes]) @@ -119,7 +118,7 @@ def to_riesz(self, poly_set): dpt_dict = self.nodes[k].deriv_dict wac_list = dpt_dict[pt] - for i in range(nbf): + for i in range(num_exp): for (w, alpha, c) in wac_list: self.matrix[k][c][i] += w*dexpansion_values[alpha][i, j] From ba29355d8a5200c6ea5c92ee025091860a7e6674 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Mon, 21 Jan 2019 10:45:59 -0600 Subject: [PATCH 25/46] make point dict and deriv dict loops look alike in to_riesz --- FIAT/dual_set.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index 276a112ee..11af03852 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -98,7 +98,7 @@ def to_riesz(self, poly_set): pt_dict = self.nodes[k].pt_dict wc_list = pt_dict[pt] - for i in range(expansion_values.shape[0]): + for i in range(num_exp]): for (w, c) in wc_list: self.matrix[k][c][i] += w*expansion_values[i, j] From d2b349b4e2f4899c16e74fd143d6c95eaa4442e6 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 1 Feb 2019 12:53:57 -0600 Subject: [PATCH 26/46] Add some documentation for to_riesz --- FIAT/dual_set.py | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index 11af03852..689c61c29 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -53,7 +53,20 @@ def get_reference_element(self): return self.ref_el def to_riesz(self, poly_set): - + """This method gives the action of the entire dual set + on each member of the expansion set underlying poly_set. + Then, applying the linear functionals of the dual set to an + arbitrary polynomial in poly_set is accomplished by matrix + multiplication.""" + + # This rather technical code queries the low-level information + # for each functional to find out where it evaluates its + # inputs and/or their derivatives. Then, it tabulates the + # expansion set one time for all the function values and + # another for all of the derivatives. This circumvents + # needing to call the to_riesz method of each functional and + # also limits the number of different calls to tabulate. + tshape = self.nodes[0].target_shape num_nodes = len(self.nodes) es = poly_set.get_expansion_set() From 41b1374a3f9b7bed6b1e0f3761cd2665c8c8432e Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 1 Feb 2019 12:59:07 -0600 Subject: [PATCH 27/46] Fix a typo and flake8 --- FIAT/dual_set.py | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index 689c61c29..d0bb8f0ae 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -58,7 +58,7 @@ def to_riesz(self, poly_set): Then, applying the linear functionals of the dual set to an arbitrary polynomial in poly_set is accomplished by matrix multiplication.""" - + # This rather technical code queries the low-level information # for each functional to find out where it evaluates its # inputs and/or their derivatives. Then, it tabulates the @@ -66,14 +66,13 @@ def to_riesz(self, poly_set): # another for all of the derivatives. This circumvents # needing to call the to_riesz method of each functional and # also limits the number of different calls to tabulate. - + tshape = self.nodes[0].target_shape num_nodes = len(self.nodes) es = poly_set.get_expansion_set() ed = poly_set.get_embedded_degree() num_exp = es.get_num_members(poly_set.get_embedded_degree()) - nbf = poly_set.get_num_members() - + riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp]) self.matrix = numpy.zeros(riesz_shape, "d") @@ -111,7 +110,7 @@ def to_riesz(self, poly_set): pt_dict = self.nodes[k].pt_dict wc_list = pt_dict[pt] - for i in range(num_exp]): + for i in range(num_exp): for (w, c) in wc_list: self.matrix[k][c][i] += w*expansion_values[i, j] From db89a0464fcd3de1b2d1e7ab4ce9da2752eb9c13 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 8 Feb 2019 12:06:44 +0000 Subject: [PATCH 28/46] Change of coordinates map implemented --- FIAT/discontinuous_serendipity.py | 61 +++++++++++++++++++++++++------ 1 file changed, 50 insertions(+), 11 deletions(-) diff --git a/FIAT/discontinuous_serendipity.py b/FIAT/discontinuous_serendipity.py index 0ec8377fc..91037072a 100644 --- a/FIAT/discontinuous_serendipity.py +++ b/FIAT/discontinuous_serendipity.py @@ -15,14 +15,22 @@ from FIAT import finite_element, polynomial_set, dual_set, functional, P0, reference_element from FIAT.discontinuous_lagrange import DiscontinuousLagrangeDualSet -from FIAT.reference_element import Point, DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron, UFCTriangle, UFCTetrahedron +from FIAT.reference_element import (Point, + DefaultLine, + UFCInterval, + UFCQuadrilateral, + UFCHexahedron, + UFCTriangle, + UFCTetrahedron, + make_affine_mapping) from FIAT.P0 import P0Dual +import numpy as np -hypercube_simplex_map = {'Point': Point(), - 'DefaultLine': DefaultLine(), - 'UFCInterval': UFCInterval(), - 'UFCQuadrilateral': UFCTriangle(), - 'UFCHexahedron': UFCTetrahedron()} +hypercube_simplex_map = {Point(): Point(), + DefaultLine(): DefaultLine(), + UFCInterval(): UFCInterval(), + UFCQuadrilateral(): UFCTriangle(), + UFCHexahedron(): UFCTetrahedron()} class DPC0Dual(P0Dual): def __init__(self, ref_el): @@ -30,7 +38,7 @@ def __init__(self, ref_el): class DPC0(finite_element.CiarletElement): def __init__(self, ref_el): - poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[type(ref_el).__name__], 0) + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) dual = DPC0Dual(ref_el) degree = 0 formdegree = ref_el.get_spatial_dimension() # n-form @@ -41,22 +49,53 @@ def __init__(self, ref_el): formdegree=formdegree) -class DiscontinuousSerendipityDualSet(DiscontinuousLagrangeDualSet): +class DiscontinuousSerendipityDualSet(dual_set.DualSet): """The dual basis for Serendipity 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): - super(DiscontinuousSerendipityDualSet, self).__init__(ref_el, degree) + entity_ids = {} + nodes = [] + + ### Change coordinates here. + v_simplex = hypercube_simplex_map[ref_el].get_vertices() + v_hypercube = ref_el.get_vertices() + v_ = list(v_hypercube[:2]) + for d in range(1, ref_el.get_dimension()): + v_.append(tuple(np.asarray(v_hypercube[d+1] + + np.average(np.asarray(v_hypercube[2**d:2**(d+1)]),axis=0)))) + A, b = make_affine_mapping(v_simplex, tuple(v_)) + + + # make nodes by getting points + # need to do this dimension-by-dimension, facet-by-facet + top = hypercube_simplex_map[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 + entity_ids[dim][entity] = [] + cur += nnodes_cur + + entity_ids[dim][0] = list(range(len(nodes))) + super(DiscontinuousSerendipityDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) class HigherOrderDiscontinuousSerendipity(finite_element.CiarletElement): """The discontinuous Serendipity finite element. It is what it is.""" def __init__(self, ref_el, degree): - poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[type(ref_el).__name__], degree) - dual = DiscontinuousSerendipityDualSet(hypercube_simplex_map[type(ref_el).__name__], degree) + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], degree) + dual = DiscontinuousSerendipityDualSet(ref_el, degree) formdegree = ref_el.get_spatial_dimension() # n-form super(HigherOrderDiscontinuousSerendipity, self).__init__(poly_set=poly_set, dual=dual, From 34eb8d3f706e8f16c7997ddfd3fa6d7a9a9024e1 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 8 Feb 2019 12:07:31 +0000 Subject: [PATCH 29/46] moved ref_el optional argument in init to the end --- FIAT/finite_element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index aba117298..15ead7129 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -120,7 +120,7 @@ class CiarletElement(FiniteElement): basis generated from polynomials encoded in a `PolynomialSet`. """ - def __init__(self, poly_set, dual, order, ref_el=None, formdegree=None, mapping="affine"): + 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) From e557bcd7f422b34592b1c414e2a04bf4593e83cd Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 8 Feb 2019 12:08:15 +0000 Subject: [PATCH 30/46] make_quadrature for hypercubes --- FIAT/quadrature.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/FIAT/quadrature.py b/FIAT/quadrature.py index 035326093..95292bb0b 100644 --- a/FIAT/quadrature.py +++ b/FIAT/quadrature.py @@ -252,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): From 000571739e3f7986269ab153664687d510d7720b Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 8 Feb 2019 12:17:03 +0000 Subject: [PATCH 31/46] Renamed discontinuous serendipity to discontinuous p --- FIAT/discontinuous_p.py | 110 ++++++++++++++++++++ test/unit/test_discontinuous_p.py | 51 +++++++++ test/unit/test_discontinuous_serendipity.py | 51 +++++++++ 3 files changed, 212 insertions(+) create mode 100644 FIAT/discontinuous_p.py create mode 100644 test/unit/test_discontinuous_p.py create mode 100644 test/unit/test_discontinuous_serendipity.py diff --git a/FIAT/discontinuous_p.py b/FIAT/discontinuous_p.py new file mode 100644 index 000000000..ee15bb815 --- /dev/null +++ b/FIAT/discontinuous_p.py @@ -0,0 +1,110 @@ +# 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 . + +from FIAT import finite_element, polynomial_set, dual_set, functional, P0, reference_element +from FIAT.reference_element import (Point, + DefaultLine, + UFCInterval, + UFCQuadrilateral, + UFCHexahedron, + UFCTriangle, + UFCTetrahedron, + make_affine_mapping) +from FIAT.P0 import P0Dual +import numpy as np + +hypercube_simplex_map = {Point(): Point(), + DefaultLine(): DefaultLine(), + UFCInterval(): UFCInterval(), + UFCQuadrilateral(): UFCTriangle(), + UFCHexahedron(): UFCTetrahedron()} + +class DPC0Dual(P0Dual): + def __init__(self, ref_el): + super(DPC0Dual, self).__init__(ref_el) + +class DPC0(finite_element.CiarletElement): + def __init__(self, ref_el): + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) + dual = DPC0Dual(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 DiscontinuousPDualSet(dual_set.DualSet): + """The dual basis for Serendipity 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. + v_simplex = hypercube_simplex_map[ref_el].get_vertices() + v_hypercube = ref_el.get_vertices() + v_ = list(v_hypercube[:2]) + for d in range(1, ref_el.get_dimension()): + v_.append(tuple(np.asarray(v_hypercube[d+1] + + np.average(np.asarray(v_hypercube[2**d:2**(d+1)]),axis=0)))) + A, b = make_affine_mapping(v_simplex, tuple(v_)) + + + # make nodes by getting points + # need to do this dimension-by-dimension, facet-by-facet + top = hypercube_simplex_map[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 + entity_ids[dim][entity] = [] + cur += nnodes_cur + + entity_ids[dim][0] = list(range(len(nodes))) + super(DiscontinuousPDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) + + +class HigherOrderDiscontinuousP(finite_element.CiarletElement): + """The discontinuous Serendipity 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 = DiscontinuousPDualSet(ref_el, degree) + formdegree = ref_el.get_spatial_dimension() # n-form + super(HigherOrderDiscontinuousP, self).__init__(poly_set=poly_set, + dual=dual, + order=degree, + ref_el=ref_el, + formdegree=formdegree) + + +def DiscontinuousP(ref_el, degree): + if degree == 0: + return DPC0(ref_el) + else: + return HigherOrderDiscontinuousP(ref_el, degree) diff --git a/test/unit/test_discontinuous_p.py b/test/unit/test_discontinuous_p.py new file mode 100644 index 000000000..8663fa44a --- /dev/null +++ b/test/unit/test_discontinuous_p.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_p import DiscontinuousP + + cell = np.array([None, 'interval', 'quadrilateral', 'hexahedron']) + s = ufc_cell(cell[dim]) + q = make_quadrature(s, degree + 1) + + fe = DiscontinuousP(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__)) diff --git a/test/unit/test_discontinuous_serendipity.py b/test/unit/test_discontinuous_serendipity.py new file mode 100644 index 000000000..3233fdd53 --- /dev/null +++ b/test/unit/test_discontinuous_serendipity.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_serendipity import DiscontinuousSerendipity + + cell = np.array([None, 'interval', 'quadrilateral', 'hexahedron']) + s = ufc_cell(cell[dim]) + q = make_quadrature(s, degree + 1) + + fe = DiscontinuousSerendipity(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__)) From c35240d5054ece42ae8a2d3b8a915e031d9156c6 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 8 Feb 2019 12:54:00 +0000 Subject: [PATCH 32/46] renamed discontinuous_p to DPC --- FIAT/discontinuous_p.py | 22 ++++---- FIAT/discontinuous_pc.py | 110 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 121 insertions(+), 11 deletions(-) create mode 100644 FIAT/discontinuous_pc.py diff --git a/FIAT/discontinuous_p.py b/FIAT/discontinuous_p.py index ee15bb815..b52231b87 100644 --- a/FIAT/discontinuous_p.py +++ b/FIAT/discontinuous_p.py @@ -48,7 +48,7 @@ def __init__(self, ref_el): formdegree=formdegree) -class DiscontinuousPDualSet(dual_set.DualSet): +class DPCDualSet(dual_set.DualSet): """The dual basis for Serendipity elements. This class works for hypercubes of any dimension. Nodes are point evaluation at equispaced points. This is the discontinuous version where @@ -86,25 +86,25 @@ def __init__(self, ref_el, degree): cur += nnodes_cur entity_ids[dim][0] = list(range(len(nodes))) - super(DiscontinuousPDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) + super(DPCDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) -class HigherOrderDiscontinuousP(finite_element.CiarletElement): +class HigherOrderDPC(finite_element.CiarletElement): """The discontinuous Serendipity 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 = DiscontinuousPDualSet(ref_el, degree) + dual = DPCDualSet(ref_el, degree) formdegree = ref_el.get_spatial_dimension() # n-form - super(HigherOrderDiscontinuousP, self).__init__(poly_set=poly_set, - dual=dual, - order=degree, - ref_el=ref_el, - formdegree=formdegree) + super(HigherOrderDPC, self).__init__(poly_set=poly_set, + dual=dual, + order=degree, + ref_el=ref_el, + formdegree=formdegree) -def DiscontinuousP(ref_el, degree): +def DPC(ref_el, degree): if degree == 0: return DPC0(ref_el) else: - return HigherOrderDiscontinuousP(ref_el, degree) + return HigherOrderDPC(ref_el, degree) diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py new file mode 100644 index 000000000..b52231b87 --- /dev/null +++ b/FIAT/discontinuous_pc.py @@ -0,0 +1,110 @@ +# 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 . + +from FIAT import finite_element, polynomial_set, dual_set, functional, P0, reference_element +from FIAT.reference_element import (Point, + DefaultLine, + UFCInterval, + UFCQuadrilateral, + UFCHexahedron, + UFCTriangle, + UFCTetrahedron, + make_affine_mapping) +from FIAT.P0 import P0Dual +import numpy as np + +hypercube_simplex_map = {Point(): Point(), + DefaultLine(): DefaultLine(), + UFCInterval(): UFCInterval(), + UFCQuadrilateral(): UFCTriangle(), + UFCHexahedron(): UFCTetrahedron()} + +class DPC0Dual(P0Dual): + def __init__(self, ref_el): + super(DPC0Dual, self).__init__(ref_el) + +class DPC0(finite_element.CiarletElement): + def __init__(self, ref_el): + poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) + dual = DPC0Dual(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 Serendipity 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. + v_simplex = hypercube_simplex_map[ref_el].get_vertices() + v_hypercube = ref_el.get_vertices() + v_ = list(v_hypercube[:2]) + for d in range(1, ref_el.get_dimension()): + v_.append(tuple(np.asarray(v_hypercube[d+1] + + np.average(np.asarray(v_hypercube[2**d:2**(d+1)]),axis=0)))) + A, b = make_affine_mapping(v_simplex, tuple(v_)) + + + # make nodes by getting points + # need to do this dimension-by-dimension, facet-by-facet + top = hypercube_simplex_map[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 + entity_ids[dim][entity] = [] + cur += nnodes_cur + + entity_ids[dim][0] = list(range(len(nodes))) + super(DPCDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) + + +class HigherOrderDPC(finite_element.CiarletElement): + """The discontinuous Serendipity 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) From 2a482e7d9220d23585fedae72a5df21b2d99a451 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Tue, 12 Feb 2019 17:44:02 +0000 Subject: [PATCH 33/46] Renamed discontinuous_pc --- FIAT/discontinuous_p.py | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/FIAT/discontinuous_p.py b/FIAT/discontinuous_p.py index b52231b87..6723311df 100644 --- a/FIAT/discontinuous_p.py +++ b/FIAT/discontinuous_p.py @@ -48,7 +48,7 @@ def __init__(self, ref_el): formdegree=formdegree) -class DPCDualSet(dual_set.DualSet): +class DPDualSet(dual_set.DualSet): """The dual basis for Serendipity elements. This class works for hypercubes of any dimension. Nodes are point evaluation at equispaced points. This is the discontinuous version where @@ -59,13 +59,20 @@ def __init__(self, ref_el, degree): 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_ = list(v_hypercube[:2]) + + # 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[d+1] + np.average(np.asarray(v_hypercube[2**d:2**(d+1)]),axis=0)))) - A, b = make_affine_mapping(v_simplex, tuple(v_)) + A, b = make_affine_mapping(v_simplex, tuple(v_)) # Make affine mapping to be used later. # make nodes by getting points @@ -86,25 +93,25 @@ def __init__(self, ref_el, degree): cur += nnodes_cur entity_ids[dim][0] = list(range(len(nodes))) - super(DPCDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) + super(DPDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) -class HigherOrderDPC(finite_element.CiarletElement): +class HigherOrderDP(finite_element.CiarletElement): """The discontinuous Serendipity 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) + dual = DPDualSet(ref_el, degree) formdegree = ref_el.get_spatial_dimension() # n-form - super(HigherOrderDPC, self).__init__(poly_set=poly_set, + super(HigherOrderDP, self).__init__(poly_set=poly_set, dual=dual, order=degree, ref_el=ref_el, formdegree=formdegree) -def DPC(ref_el, degree): +def DP(ref_el, degree): if degree == 0: return DPC0(ref_el) else: - return HigherOrderDPC(ref_el, degree) + return HigherOrderDP(ref_el, degree) From d538802e3e89c6470394f29457f759134f12d0d6 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Tue, 12 Feb 2019 17:44:57 +0000 Subject: [PATCH 34/46] Renamed discontinuous_pc --- test/unit/test_discontinuous_p.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/test_discontinuous_p.py b/test/unit/test_discontinuous_p.py index 8663fa44a..5b15ddfe0 100644 --- a/test/unit/test_discontinuous_p.py +++ b/test/unit/test_discontinuous_p.py @@ -29,13 +29,13 @@ 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_p import DiscontinuousP + from FIAT.discontinuous_p import DP cell = np.array([None, 'interval', 'quadrilateral', 'hexahedron']) s = ufc_cell(cell[dim]) q = make_quadrature(s, degree + 1) - fe = DiscontinuousP(s, degree) + fe = DP(s, degree) tab = fe.tabulate(0, q.pts)[(0,) * dim] for test_degree in range(degree + 1): From 73fda42f954d5a43930d3c92e010a1fcf89c9e5b Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Thu, 14 Feb 2019 10:15:35 +0000 Subject: [PATCH 35/46] Renamed test file for DPC --- test/unit/test_discontinuous_pc.py | 51 ++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 test/unit/test_discontinuous_pc.py 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__)) From 5204a813d038d7f0720dc673980abfb979ed64b9 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Thu, 14 Feb 2019 10:16:55 +0000 Subject: [PATCH 36/46] Added DPC to init file --- FIAT/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 8e60bb270..622ab65d1 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -15,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 From c13e4ced06c3293293f6c710c499ed6c179f9e85 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Thu, 14 Feb 2019 10:17:28 +0000 Subject: [PATCH 37/46] Renamed DP to DPC --- FIAT/discontinuous_pc.py | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index b52231b87..552c5968a 100644 --- a/FIAT/discontinuous_pc.py +++ b/FIAT/discontinuous_pc.py @@ -59,18 +59,26 @@ def __init__(self, ref_el, degree): 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_ = list(v_hypercube[:2]) + + # 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[d+1] + np.average(np.asarray(v_hypercube[2**d:2**(d+1)]),axis=0)))) - A, b = make_affine_mapping(v_simplex, tuple(v_)) + 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): @@ -82,8 +90,9 @@ def __init__(self, ref_el, degree): for x in pts_cur] nnodes_cur = len(nodes_cur) nodes += nodes_cur - entity_ids[dim][entity] = [] 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, hypercube_simplex_map[ref_el], entity_ids) From fc0de5f9be61cdabf043012912debcbe60166dac Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Thu, 14 Feb 2019 10:18:30 +0000 Subject: [PATCH 38/46] Deleted old file --- FIAT/discontinuous_p.py | 117 ---------------------------------------- 1 file changed, 117 deletions(-) delete mode 100644 FIAT/discontinuous_p.py diff --git a/FIAT/discontinuous_p.py b/FIAT/discontinuous_p.py deleted file mode 100644 index 6723311df..000000000 --- a/FIAT/discontinuous_p.py +++ /dev/null @@ -1,117 +0,0 @@ -# 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 . - -from FIAT import finite_element, polynomial_set, dual_set, functional, P0, reference_element -from FIAT.reference_element import (Point, - DefaultLine, - UFCInterval, - UFCQuadrilateral, - UFCHexahedron, - UFCTriangle, - UFCTetrahedron, - make_affine_mapping) -from FIAT.P0 import P0Dual -import numpy as np - -hypercube_simplex_map = {Point(): Point(), - DefaultLine(): DefaultLine(), - UFCInterval(): UFCInterval(), - UFCQuadrilateral(): UFCTriangle(), - UFCHexahedron(): UFCTetrahedron()} - -class DPC0Dual(P0Dual): - def __init__(self, ref_el): - super(DPC0Dual, self).__init__(ref_el) - -class DPC0(finite_element.CiarletElement): - def __init__(self, ref_el): - poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) - dual = DPC0Dual(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 DPDualSet(dual_set.DualSet): - """The dual basis for Serendipity 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_ = list(v_hypercube[:2]) - - # 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[d+1] - + np.average(np.asarray(v_hypercube[2**d:2**(d+1)]),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() - - 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 - entity_ids[dim][entity] = [] - cur += nnodes_cur - - entity_ids[dim][0] = list(range(len(nodes))) - super(DPDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) - - -class HigherOrderDP(finite_element.CiarletElement): - """The discontinuous Serendipity 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 = DPDualSet(ref_el, degree) - formdegree = ref_el.get_spatial_dimension() # n-form - super(HigherOrderDP, self).__init__(poly_set=poly_set, - dual=dual, - order=degree, - ref_el=ref_el, - formdegree=formdegree) - - -def DP(ref_el, degree): - if degree == 0: - return DPC0(ref_el) - else: - return HigherOrderDP(ref_el, degree) From f1ddb41343af2dc23075de0c6d51d40679c292c4 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Thu, 14 Feb 2019 10:18:46 +0000 Subject: [PATCH 39/46] Deleted old file --- test/unit/test_discontinuous_p.py | 51 ------------------------------- 1 file changed, 51 deletions(-) delete mode 100644 test/unit/test_discontinuous_p.py diff --git a/test/unit/test_discontinuous_p.py b/test/unit/test_discontinuous_p.py deleted file mode 100644 index 5b15ddfe0..000000000 --- a/test/unit/test_discontinuous_p.py +++ /dev/null @@ -1,51 +0,0 @@ -# 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_p import DP - - cell = np.array([None, 'interval', 'quadrilateral', 'hexahedron']) - s = ufc_cell(cell[dim]) - q = make_quadrature(s, degree + 1) - - fe = DP(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__)) From e937846680d0898cba21908491eb65e389989727 Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Fri, 15 Feb 2019 13:58:27 +0000 Subject: [PATCH 40/46] Address review comments --- FIAT/dual_set.py | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index d0bb8f0ae..18884b582 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -75,7 +75,7 @@ def to_riesz(self, poly_set): riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp]) - self.matrix = numpy.zeros(riesz_shape, "d") + result = numpy.zeros(riesz_shape, "d") # let's amalgamate the pt_dict and deriv_dicts of all the # functionals so we can tabulate the basis functions twice only @@ -87,17 +87,11 @@ def to_riesz(self, poly_set): for i, ell in enumerate(self.nodes): for pt in ell.pt_dict: - if pt in pts_to_ells: - pts_to_ells[pt].append(i) - else: - pts_to_ells[pt] = [i] + pts_to_ells.setdefault(pt, []).append(i) for i, ell in enumerate(self.nodes): for pt in ell.deriv_dict: - if pt in dpts_to_ells: - dpts_to_ells[pt].append(i) - else: - dpts_to_ells[pt] = [i] + dpts_to_ells.setdefault(pt, []).append(i) # Now tabulate pts = list(pts_to_ells.keys()) @@ -112,7 +106,7 @@ def to_riesz(self, poly_set): for i in range(num_exp): for (w, c) in wc_list: - self.matrix[k][c][i] += w*expansion_values[i, j] + result[k][c][i] += w*expansion_values[i, j] max_deriv_order = max([ell.max_deriv_order for ell in self.nodes]) if max_deriv_order > 0: @@ -132,6 +126,6 @@ def to_riesz(self, poly_set): for i in range(num_exp): for (w, alpha, c) in wac_list: - self.matrix[k][c][i] += w*dexpansion_values[alpha][i, j] + result[k][c][i] += w*dexpansion_values[alpha][i, j] - return self.matrix + return result From 358a4cfde8bc597b9e20597cae0be0a7bbd64dc5 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Sat, 16 Feb 2019 19:22:03 +0000 Subject: [PATCH 41/46] Removed unused files --- FIAT/discontinuous_serendipity.py | 111 -------------------- test/unit/test_discontinuous_serendipity.py | 51 --------- 2 files changed, 162 deletions(-) delete mode 100644 FIAT/discontinuous_serendipity.py delete mode 100644 test/unit/test_discontinuous_serendipity.py diff --git a/FIAT/discontinuous_serendipity.py b/FIAT/discontinuous_serendipity.py deleted file mode 100644 index 91037072a..000000000 --- a/FIAT/discontinuous_serendipity.py +++ /dev/null @@ -1,111 +0,0 @@ -# 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 . - -from FIAT import finite_element, polynomial_set, dual_set, functional, P0, reference_element -from FIAT.discontinuous_lagrange import DiscontinuousLagrangeDualSet -from FIAT.reference_element import (Point, - DefaultLine, - UFCInterval, - UFCQuadrilateral, - UFCHexahedron, - UFCTriangle, - UFCTetrahedron, - make_affine_mapping) -from FIAT.P0 import P0Dual -import numpy as np - -hypercube_simplex_map = {Point(): Point(), - DefaultLine(): DefaultLine(), - UFCInterval(): UFCInterval(), - UFCQuadrilateral(): UFCTriangle(), - UFCHexahedron(): UFCTetrahedron()} - -class DPC0Dual(P0Dual): - def __init__(self, ref_el): - super(DPC0Dual, self).__init__(ref_el) - -class DPC0(finite_element.CiarletElement): - def __init__(self, ref_el): - poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) - dual = DPC0Dual(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 DiscontinuousSerendipityDualSet(dual_set.DualSet): - """The dual basis for Serendipity 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. - v_simplex = hypercube_simplex_map[ref_el].get_vertices() - v_hypercube = ref_el.get_vertices() - v_ = list(v_hypercube[:2]) - for d in range(1, ref_el.get_dimension()): - v_.append(tuple(np.asarray(v_hypercube[d+1] - + np.average(np.asarray(v_hypercube[2**d:2**(d+1)]),axis=0)))) - A, b = make_affine_mapping(v_simplex, tuple(v_)) - - - # make nodes by getting points - # need to do this dimension-by-dimension, facet-by-facet - top = hypercube_simplex_map[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 - entity_ids[dim][entity] = [] - cur += nnodes_cur - - entity_ids[dim][0] = list(range(len(nodes))) - super(DiscontinuousSerendipityDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) - - -class HigherOrderDiscontinuousSerendipity(finite_element.CiarletElement): - """The discontinuous Serendipity 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 = DiscontinuousSerendipityDualSet(ref_el, degree) - formdegree = ref_el.get_spatial_dimension() # n-form - super(HigherOrderDiscontinuousSerendipity, self).__init__(poly_set=poly_set, - dual=dual, - order=degree, - ref_el=ref_el, - formdegree=formdegree) - - -def DiscontinuousSerendipity(ref_el, degree): - if degree == 0: - return DPC0(ref_el) - else: - return HigherOrderDiscontinuousSerendipity(ref_el, degree) diff --git a/test/unit/test_discontinuous_serendipity.py b/test/unit/test_discontinuous_serendipity.py deleted file mode 100644 index 3233fdd53..000000000 --- a/test/unit/test_discontinuous_serendipity.py +++ /dev/null @@ -1,51 +0,0 @@ -# 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_serendipity import DiscontinuousSerendipity - - cell = np.array([None, 'interval', 'quadrilateral', 'hexahedron']) - s = ufc_cell(cell[dim]) - q = make_quadrature(s, degree + 1) - - fe = DiscontinuousSerendipity(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__)) From 974c20980561311994322d2551f92794602aa539 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Sat, 16 Feb 2019 19:29:17 +0000 Subject: [PATCH 42/46] Edited docstrings --- FIAT/discontinuous_pc.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index 552c5968a..f092662c8 100644 --- a/FIAT/discontinuous_pc.py +++ b/FIAT/discontinuous_pc.py @@ -49,7 +49,7 @@ def __init__(self, ref_el): class DPCDualSet(dual_set.DualSet): - """The dual basis for Serendipity elements. This class works for + """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""" @@ -99,7 +99,7 @@ def __init__(self, ref_el, degree): class HigherOrderDPC(finite_element.CiarletElement): - """The discontinuous Serendipity finite element. It is what it is.""" + """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) From 54c4a4ad932c9beae4110250f66aae0fab6d1e47 Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Tue, 19 Feb 2019 16:47:38 +0000 Subject: [PATCH 43/46] Revert "Merged in branch 'toriesz' (pull request #54)" This reverts commit 1375329334c09416148c5f23134cf5b8e7fdaada, reversing changes made to daacf42c19a83843ffcde3dc18fc0d5922feebcb. This breaks the Bell element, and also FacetBubble. --- FIAT/dual_set.py | 73 ++--------------------- FIAT/functional.py | 90 ++++++++++++++--------------- bitbucket-pipelines.yml | 4 +- test/regression/test_regression.py | 2 +- test/unit/test_fiat.py | 64 ++++++++++---------- test/unit/test_quadrature.py | 40 ++++--------- test/unit/test_reference_element.py | 18 +++--- 7 files changed, 102 insertions(+), 189 deletions(-) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index 18884b582..b9e2e74e2 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -16,9 +16,6 @@ # along with FIAT. If not, see . import numpy -import collections - -from FIAT import polynomial_set class DualSet(object): @@ -53,79 +50,17 @@ def get_reference_element(self): return self.ref_el def to_riesz(self, poly_set): - """This method gives the action of the entire dual set - on each member of the expansion set underlying poly_set. - Then, applying the linear functionals of the dual set to an - arbitrary polynomial in poly_set is accomplished by matrix - multiplication.""" - - # This rather technical code queries the low-level information - # for each functional to find out where it evaluates its - # inputs and/or their derivatives. Then, it tabulates the - # expansion set one time for all the function values and - # another for all of the derivatives. This circumvents - # needing to call the to_riesz method of each functional and - # also limits the number of different calls to tabulate. tshape = self.nodes[0].target_shape num_nodes = len(self.nodes) es = poly_set.get_expansion_set() - ed = poly_set.get_embedded_degree() num_exp = es.get_num_members(poly_set.get_embedded_degree()) riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp]) - result = numpy.zeros(riesz_shape, "d") - - # let's amalgamate the pt_dict and deriv_dicts of all the - # functionals so we can tabulate the basis functions twice only - # (once on pts and once on derivatives) - - # Need: dictionary mapping pts to which functionals they come from - pts_to_ells = collections.OrderedDict() - dpts_to_ells = collections.OrderedDict() - - for i, ell in enumerate(self.nodes): - for pt in ell.pt_dict: - pts_to_ells.setdefault(pt, []).append(i) - - for i, ell in enumerate(self.nodes): - for pt in ell.deriv_dict: - dpts_to_ells.setdefault(pt, []).append(i) - - # Now tabulate - pts = list(pts_to_ells.keys()) - expansion_values = es.tabulate(ed, pts) - - for j, pt in enumerate(pts): - which_ells = pts_to_ells[pt] - - for k in which_ells: - pt_dict = self.nodes[k].pt_dict - wc_list = pt_dict[pt] - - for i in range(num_exp): - for (w, c) in wc_list: - result[k][c][i] += w*expansion_values[i, j] - - max_deriv_order = max([ell.max_deriv_order for ell in self.nodes]) - if max_deriv_order > 0: - dpts = list(dpts_to_ells.keys()) - # It's easiest/most efficient to get derivatives of the - # expansion set through the polynomial set interface. - # This is creating a short-lived set to do just this. - expansion = polynomial_set.ONPolynomialSet(self.ref_el, ed) - dexpansion_values = expansion.tabulate(dpts, max_deriv_order) - - for j, pt in enumerate(dpts): - which_ells = dpts_to_ells[pt] - - for k in which_ells: - dpt_dict = self.nodes[k].deriv_dict - wac_list = dpt_dict[pt] + self.mat = numpy.zeros(riesz_shape, "d") - for i in range(num_exp): - for (w, alpha, c) in wac_list: - result[k][c][i] += w*dexpansion_values[alpha][i, j] + for i in range(len(self.nodes)): + self.mat[i][:] = self.nodes[i].to_riesz(poly_set) - return result + return self.mat diff --git a/FIAT/functional.py b/FIAT/functional.py index 6eb527a5e..8d66327a2 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -26,8 +26,6 @@ import numpy import sympy -from FIAT import polynomial_set - def index_iterator(shp): """Constructs a generator iterating over all indices in @@ -44,6 +42,11 @@ def index_iterator(shp): for foo in index_iterator(shp_foo): yield [i] + foo +# also put in a "jet_dict" that maps +# pt --> {wt, multiindex, comp} +# the multiindex is an iterable of nonnegative +# integers + class Functional(object): """Class implementing an abstract functional. @@ -91,6 +94,7 @@ def get_type_tag(self): normal component, which is probably handy for clients of FIAT""" return self.functional_type + # overload me in subclasses to make life easier!! def to_riesz(self, poly_set): """Constructs an array representation of the functional over the base of the given polynomial_set so that f(phi) for any @@ -99,20 +103,17 @@ def to_riesz(self, poly_set): ed = poly_set.get_embedded_degree() pt_dict = self.get_point_dict() - nbf = poly_set.get_num_members() - pts = list(pt_dict.keys()) # bfs is matrix that is pdim rows by num_pts cols # where pdim is the polynomial dimension bfs = es.tabulate(ed, pts) - npts = len(pts) result = numpy.zeros(poly_set.coeffs.shape[1:], "d") # loop over points - for j in range(npts): + for j in range(len(pts)): pt_cur = pts[j] wc_list = pt_dict[pt_cur] @@ -122,23 +123,7 @@ def to_riesz(self, poly_set): result[c][i] += w * bfs[i, j] if self.deriv_dict: - dpt_dict = self.deriv_dict - - # this makes things quicker since it uses dmats after - # instantiation - es_foo = polynomial_set.ONPolynomialSet(self.ref_el, ed) - dpts = list(dpt_dict.keys()) - - # figure out the maximum order of derivative being taken - dbfs = es_foo.tabulate(dpts, self.max_deriv_order) - - ndpts = len(dpts) - for j in range(ndpts): - dpt_cur = dpts[j] - wac_list = dpt_dict[dpt_cur] - for i in range(nbf): - for (w, alpha, c) in wac_list: - result[c][i] += w * dbfs[alpha][i, j] + raise NotImplementedError("Generic to_riesz implementation does not support derivatives") return result @@ -187,8 +172,8 @@ class PointDerivative(Functional): functions at a particular point x.""" def __init__(self, ref_el, x, alpha): - dpt_dict = {x: [(1.0, tuple(alpha), tuple())]} - self.alpha = tuple(alpha) + dpt_dict = {x: [(1.0, alpha, tuple())]} + self.alpha = alpha self.order = sum(self.alpha) Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointDeriv") @@ -206,6 +191,25 @@ def __call__(self, fn): return sympy.diff(fn(X), *dvars).evalf(subs=dict(zip(dX, x))) + def to_riesz(self, poly_set): + x = list(self.deriv_dict.keys())[0] + + X = sympy.DeferredVector('x') + dx = numpy.asarray([X[i] for i in range(len(x))]) + + es = poly_set.get_expansion_set() + ed = poly_set.get_embedded_degree() + + bfs = es.tabulate(ed, [dx])[:, 0] + + # Expand the multi-index as a series of variables to + # differentiate with respect to. + dvars = tuple(d for d, a in zip(dx, self.alpha) + for count in range(a)) + + return numpy.asarray([sympy.lambdify(X, sympy.diff(b, *dvars))(x) + for b in bfs]) + class PointNormalDerivative(Functional): @@ -219,35 +223,27 @@ def __init__(self, ref_el, facet_no, pt): alpha = [0] * sd alpha[i] = 1 alphas.append(alpha) - dpt_dict = {pt: [(n[i], tuple(alphas[i]), tuple()) for i in range(sd)]} + dpt_dict = {pt: [(n[i], alphas[i], tuple()) for i in range(sd)]} Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv") + def to_riesz(self, poly_set): + x = list(self.deriv_dict.keys())[0] -class PointNormalSecondDerivative(Functional): + X = sympy.DeferredVector('x') + dx = numpy.asarray([X[i] for i in range(len(x))]) - def __init__(self, ref_el, facet_no, pt): - n = ref_el.compute_normal(facet_no) - self.n = n - sd = ref_el.get_spatial_dimension() - tau = numpy.zeros((sd*(sd+1)//2,)) + es = poly_set.get_expansion_set() + ed = poly_set.get_embedded_degree() - alphas = [] - cur = 0 - for i in range(sd): - for j in range(i, sd): - alpha = [0] * sd - alpha[i] += 1 - alpha[j] += 1 - alphas.append(tuple(alpha)) - tau[cur] = n[i]*n[j] - cur += 1 - - self.tau = tau - self.alphas = alphas - dpt_dict = {pt: [(n[i], alphas[i], tuple()) for i in range(sd)]} + bfs = es.tabulate(ed, [dx])[:, 0] - Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv") + # We need the gradient dotted with the normal. + return numpy.asarray( + [sympy.lambdify( + X, sum([sympy.diff(b, dxi)*ni + for dxi, ni in zip(dx, self.n)]))(x) + for b in bfs]) class IntegralMoment(Functional): diff --git a/bitbucket-pipelines.yml b/bitbucket-pipelines.yml index e3488e10e..a32b26f1c 100644 --- a/bitbucket-pipelines.yml +++ b/bitbucket-pipelines.yml @@ -6,9 +6,9 @@ pipelines: script: - apt-get update - apt-get install -y - git python3-minimal python3-pip python3-setuptools python3-wheel + git python3-minimal python3-pip python3-pytest python3-setuptools python3-wheel --no-install-recommends - - pip3 install flake8 pytest --upgrade + - pip3 install flake8 - pip3 install . - python3 -m flake8 - DATA_REPO_GIT="" python3 -m pytest -v test/ diff --git a/test/regression/test_regression.py b/test/regression/test_regression.py index 4d8afc572..f2f510f41 100644 --- a/test/regression/test_regression.py +++ b/test/regression/test_regression.py @@ -277,7 +277,7 @@ def _perform_test(family, dim, degree, reference_table): reference = json.load(open(filename, "r"), object_hook=json_numpy_obj_hook) for test_case in test_cases: family, dim, degree = test_case - _perform_test(family, dim, degree, reference[str(test_case)]) + yield _perform_test, family, dim, degree, reference[str(test_case)] # Update references if missing except (IOError, KeyError) as e: diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index e949d0468..eec0bd817 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -203,42 +203,42 @@ def __init__(self, a, b): # MixedElement made of nodal elements should be nodal, but its API # is currently just broken. - pytest.param("MixedElement([" - " DiscontinuousLagrange(T, 1)," - " RaviartThomas(T, 2)" - "])", marks=xfail_impl), + xfail_impl("MixedElement([" + " DiscontinuousLagrange(T, 1)," + " RaviartThomas(T, 2)" + "])"), # Following element do not bother implementing get_nodal_basis # so the test would need to be rewritten using tabulate - pytest.param("TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2))", marks=xfail_impl), - pytest.param("Hdiv(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))", marks=xfail_impl), - pytest.param("Hcurl(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))", marks=xfail_impl), - pytest.param("HDivTrace(T, 1)", marks=xfail_impl), - pytest.param("EnrichedElement(" - "Hdiv(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " - "Hdiv(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" - ")", marks=xfail_impl), - pytest.param("EnrichedElement(" - "Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " - "Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" - ")", marks=xfail_impl), + xfail_impl("TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2))"), + xfail_impl("Hdiv(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"), + xfail_impl("Hcurl(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"), + xfail_impl("HDivTrace(T, 1)"), + xfail_impl("EnrichedElement(" + "Hdiv(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " + "Hdiv(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" + ")"), + xfail_impl("EnrichedElement(" + "Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), " + "Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))" + ")"), # Following elements are checked using tabulate - pytest.param("HDivTrace(T, 0)", marks=xfail_impl), - pytest.param("HDivTrace(T, 1)", marks=xfail_impl), - pytest.param("HDivTrace(T, 2)", marks=xfail_impl), - pytest.param("HDivTrace(T, 3)", marks=xfail_impl), - pytest.param("HDivTrace(S, 0)", marks=xfail_impl), - pytest.param("HDivTrace(S, 1)", marks=xfail_impl), - pytest.param("HDivTrace(S, 2)", marks=xfail_impl), - pytest.param("HDivTrace(S, 3)", marks=xfail_impl), - pytest.param("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))", marks=xfail_impl), - pytest.param("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))", marks=xfail_impl), - pytest.param("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))", marks=xfail_impl), - pytest.param("TensorProductElement(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)), Lagrange(I, 2))", marks=xfail_impl), - pytest.param("FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)))", marks=xfail_impl), - pytest.param("FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)))", marks=xfail_impl), - pytest.param("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))", marks=xfail_impl), - pytest.param("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))", marks=xfail_impl), + xfail_impl("HDivTrace(T, 0)"), + xfail_impl("HDivTrace(T, 1)"), + xfail_impl("HDivTrace(T, 2)"), + xfail_impl("HDivTrace(T, 3)"), + xfail_impl("HDivTrace(S, 0)"), + xfail_impl("HDivTrace(S, 1)"), + xfail_impl("HDivTrace(S, 2)"), + xfail_impl("HDivTrace(S, 3)"), + xfail_impl("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))"), + xfail_impl("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))"), + xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))"), + xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)), Lagrange(I, 2))"), + xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)))"), + xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)))"), + xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))"), + xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))"), ] diff --git a/test/unit/test_quadrature.py b/test/unit/test_quadrature.py index fdcf88e66..68c5c7e1c 100644 --- a/test/unit/test_quadrature.py +++ b/test/unit/test_quadrature.py @@ -134,39 +134,21 @@ def test_create_quadrature_extr_quadrilateral(extr_quadrilateral, basedeg, extrd (2**(basedeg + 2) - 2) / ((basedeg + 1)*(basedeg + 2)) * 1/(extrdeg + 1)) -def test_invalid_quadrature_degree_interval(interval, scheme): +@pytest.mark.parametrize("cell", [interval(), + triangle(), + tetrahedron(), + quadrilateral()]) +def test_invalid_quadrature_degree(cell, scheme): with pytest.raises(ValueError): - FIAT.create_quadrature(interval, -1, scheme) + FIAT.create_quadrature(cell, -1, scheme) -def test_invalid_quadrature_degree_tri(triangle, scheme): +@pytest.mark.parametrize("cell", [extr_interval(), + extr_triangle(), + extr_quadrilateral()]) +def test_invalid_quadrature_degree_tensor_prod(cell): with pytest.raises(ValueError): - FIAT.create_quadrature(triangle, -1, scheme) - - -def test_invalid_quadrature_degree_quad(quadrilateral, scheme): - with pytest.raises(ValueError): - FIAT.create_quadrature(quadrilateral, -1, scheme) - - -def test_invalid_quadrature_degree_tet(tetrahedron, scheme): - with pytest.raises(ValueError): - FIAT.create_quadrature(tetrahedron, -1, scheme) - - -def test_invalid_quadrature_degree_tensor_prod_interval(extr_interval): - with pytest.raises(ValueError): - FIAT.create_quadrature(extr_interval, (-1, -1)) - - -def test_invalid_quadrature_degree_tensor_prod_tri(extr_triangle): - with pytest.raises(ValueError): - FIAT.create_quadrature(extr_triangle, (-1, -1)) - - -def test_invalid_quadrature_degree_tensor_prod_quad(extr_quadrilateral): - with pytest.raises(ValueError): - FIAT.create_quadrature(extr_quadrilateral, (-1, -1)) + FIAT.create_quadrature(cell, (-1, -1)) def test_tensor_product_composition(interval, triangle, extr_triangle, scheme): diff --git a/test/unit/test_reference_element.py b/test/unit/test_reference_element.py index a509b6c64..edc1fea65 100644 --- a/test/unit/test_reference_element.py +++ b/test/unit/test_reference_element.py @@ -39,8 +39,8 @@ @pytest.mark.parametrize(('cell', 'connectivity'), [(tetrahedron, ufc_tetrahedron_21connectivity), (hexahedron, ufc_hexahedron_21connectivity), - pytest.param(triangle_x_interval, [], marks=pytest.mark.xfail), - pytest.param(quadrilateral_x_interval, [], marks=pytest.mark.xfail)]) + pytest.mark.xfail((triangle_x_interval, [])), + pytest.mark.xfail((quadrilateral_x_interval, []))]) def test_ufc_connectivity_21(cell, connectivity): """Check face-edge connectivity builds what UFC expects. This is only non-trivial case ; the rest is x-0 and D-x, @@ -51,9 +51,9 @@ def test_ufc_connectivity_21(cell, connectivity): @pytest.mark.parametrize('cell', [point, interval, triangle, tetrahedron, quadrilateral, hexahedron, - pytest.param(interval_x_interval, marks=pytest.mark.xfail), - pytest.param(triangle_x_interval, marks=pytest.mark.xfail), - pytest.param(quadrilateral_x_interval, marks=pytest.mark.xfail)]) + pytest.mark.xfail(interval_x_interval), + pytest.mark.xfail(triangle_x_interval), + pytest.mark.xfail(quadrilateral_x_interval)]) def test_ufc_connectivity_x0(cell): """Check x-0 connectivity is just what get_topology gives""" for dim0 in range(cell.get_spatial_dimension()+1): @@ -66,9 +66,9 @@ def test_ufc_connectivity_x0(cell): @pytest.mark.parametrize('cell', [point, interval, triangle, tetrahedron, quadrilateral, hexahedron, - pytest.param(interval_x_interval, marks=pytest.mark.xfail), - pytest.param(triangle_x_interval, marks=pytest.mark.xfail), - pytest.param(quadrilateral_x_interval, marks=pytest.mark.xfail)]) + pytest.mark.xfail(interval_x_interval), + pytest.mark.xfail(triangle_x_interval), + pytest.mark.xfail(quadrilateral_x_interval)]) def test_ufc_connectivity_Dx(cell): """Check D-x connectivity is just [(0,1,2,...)]""" D = cell.get_spatial_dimension() @@ -79,7 +79,7 @@ def test_ufc_connectivity_Dx(cell): @pytest.mark.parametrize(('cell', 'volume'), - [pytest.param(point, 1, marks=pytest.mark.xfail), + [pytest.mark.xfail((point, 1)), (interval, 1), (triangle, 1/2), (quadrilateral, 1), From 640fb21e8141c5c187f93403f1920da93542887e Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Tue, 19 Feb 2019 18:29:54 +0000 Subject: [PATCH 44/46] flake8 compliant --- FIAT/__init__.py | 1 + FIAT/discontinuous_pc.py | 13 +++++++------ 2 files changed, 8 insertions(+), 6 deletions(-) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 622ab65d1..cc399210a 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -56,6 +56,7 @@ "FacetBubble": FacetBubble, "Crouzeix-Raviart": CrouzeixRaviart, "Discontinuous Lagrange": DiscontinuousLagrange, + "DPC": DPC, "Discontinuous Taylor": DiscontinuousTaylor, "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas, "Hermite": CubicHermite, diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index f092662c8..2450fe4cf 100644 --- a/FIAT/discontinuous_pc.py +++ b/FIAT/discontinuous_pc.py @@ -13,7 +13,7 @@ # You should have received a copy of the GNU Lesser General Public License # along with FIAT. If not, see . -from FIAT import finite_element, polynomial_set, dual_set, functional, P0, reference_element +from FIAT import finite_element, polynomial_set, dual_set, functional from FIAT.reference_element import (Point, DefaultLine, UFCInterval, @@ -31,10 +31,12 @@ UFCQuadrilateral(): UFCTriangle(), UFCHexahedron(): UFCTetrahedron()} + class DPC0Dual(P0Dual): def __init__(self, ref_el): super(DPC0Dual, self).__init__(ref_el) + class DPC0(finite_element.CiarletElement): def __init__(self, ref_el): poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) @@ -58,7 +60,7 @@ def __init__(self, ref_el, degree): entity_ids = {} nodes = [] - ### Change coordinates here. + # 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. @@ -70,10 +72,9 @@ def __init__(self, ref_el, degree): # 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[d+1] - + np.average(np.asarray(v_hypercube[2**d:2**(d+1)]),axis=0)))) - A, b = make_affine_mapping(v_simplex, tuple(v_)) # Make affine mapping to be used later. - + v_.append(tuple(np.asarray(v_hypercube[d+1] + + np.average(np.asarray(v_hypercube[2**d:2**(d+1)]), 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 From 290dbafe65953f23b87ce1806114427c2c3d2cfd Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Wed, 20 Feb 2019 15:32:10 +0000 Subject: [PATCH 45/46] Deleted DPC0Dual --- FIAT/discontinuous_pc.py | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index 2450fe4cf..2811ced91 100644 --- a/FIAT/discontinuous_pc.py +++ b/FIAT/discontinuous_pc.py @@ -32,15 +32,10 @@ UFCHexahedron(): UFCTetrahedron()} -class DPC0Dual(P0Dual): - def __init__(self, ref_el): - super(DPC0Dual, self).__init__(ref_el) - - class DPC0(finite_element.CiarletElement): def __init__(self, ref_el): poly_set = polynomial_set.ONPolynomialSet(hypercube_simplex_map[ref_el], 0) - dual = DPC0Dual(ref_el) + dual = P0Dual(ref_el) degree = 0 formdegree = ref_el.get_spatial_dimension() # n-form super(DPC0, self).__init__(poly_set=poly_set, From 67c1cd0546be430bf001aee54e4a8b94820a3071 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 5 Apr 2019 12:11:03 +0100 Subject: [PATCH 46/46] fixed ref_el and coordinate change --- FIAT/discontinuous_pc.py | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/FIAT/discontinuous_pc.py b/FIAT/discontinuous_pc.py index 2811ced91..8c5683109 100644 --- a/FIAT/discontinuous_pc.py +++ b/FIAT/discontinuous_pc.py @@ -1,3 +1,5 @@ +# 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 @@ -12,8 +14,11 @@ # # 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, @@ -22,7 +27,6 @@ UFCTriangle, UFCTetrahedron, make_affine_mapping) -from FIAT.P0 import P0Dual import numpy as np hypercube_simplex_map = {Point(): Point(), @@ -61,14 +65,14 @@ def __init__(self, ref_el, degree): # Vertices of the reference element. v_hypercube = ref_el.get_vertices() # For the mapping, first two vertices are unchanged in all dimensions. - v_ = list(v_hypercube[:2]) + 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[d+1] + - np.average(np.asarray(v_hypercube[2**d:2**(d+1)]), axis=0)))) + 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 @@ -91,7 +95,7 @@ def __init__(self, ref_el, degree): entity_ids[dim][entity] = [] entity_ids[dim][0] = list(range(len(nodes))) - super(DPCDualSet, self).__init__(nodes, hypercube_simplex_map[ref_el], entity_ids) + super(DPCDualSet, self).__init__(nodes, ref_el, entity_ids) class HigherOrderDPC(finite_element.CiarletElement):