From 1d8589e5aa0fc932d1c56792aed2defe98081528 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Tue, 25 Jun 2019 23:30:39 +0100 Subject: [PATCH 01/23] started bdmc element --- FIAT/brezzi_douglas_marini_cube.py | 65 ++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 FIAT/brezzi_douglas_marini_cube.py diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py new file mode 100644 index 000000000..89a7796b6 --- /dev/null +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -0,0 +1,65 @@ +# Copyright (C) 2019 Cyrus Cheng (Imperial College London) +# +# This file is part of FIAT. +# +# FIAT is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# FIAT is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with FIAT. If not, see . +# +# Modified by David A. Ham (david.ham@imperial.ac.uk), 2019 + +from sympy import symbols, legendre, Array, diff +import numpy as np +from FIAT.finite_element import FiniteElement +from FIAT.lagrange import Lagrange +from FIAT.dual_set import make_entity_closure_ids +from FIAT.polynomial_set import mis +from FIAT.serendipity import tr + +x, y, z = symbols('x y z') +variables = (x, y, z) +leg = legendre + + +class BrezziDouglasMariniCubeEdge(FiniteElement): + def __init__(self, ref_el, degree): + if degree < 1: + raise Exception("BDMce_k elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + if dim != 2: + raise Exception("BDMce_k elements only valid for dimension 2") + + flat_topology = flat_el.get_topology() + + verts = flat_el.get_vertices() + + dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) + dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) + x_mid = 2*x-(verts[-1][0] + verts[0][0]) + y_mid = 2*y-(verts[-1][1] + verts[0][1]) + + EL = e_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + ETL = e_tilda_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + + +def e_lambda_1_2d(i, dx, dy, x_mid, y_mid): + EL = tuple([(0, y_mid**j*a) for a in dx for j in range(i)] + [(x_mid**j*a, 0) for a in dy for j in range(i)]) + + return EL + +def e_tilda_lambda_1_2d(r, dx, dy, x_mid, y_mid): + ETL = tuple([(y_mid**(r-1)*dy[0]*dy[1], (r+1)*y_mid**r*a) for a in dx] + + [((r+1)*x_mid**r*a, x_mid**(r-1)*dx[0]*dx[1]) for a in dy]) + + return ETL From 149247509373f7d4ba5a927e9e2e4d02f01c285e Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Tue, 25 Jun 2019 23:56:50 +0100 Subject: [PATCH 02/23] started implementing entity_dofs --- FIAT/brezzi_douglas_marini_cube.py | 40 ++++++++++++++++++++++++------ 1 file changed, 33 insertions(+), 7 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 89a7796b6..ce6cd9d85 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -50,16 +50,42 @@ def __init__(self, ref_el, degree): y_mid = 2*y-(verts[-1][1] + verts[0][1]) EL = e_lambda_1_2d(degree, dx, dy, x_mid, y_mid) - ETL = e_tilda_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + FL = f_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + bdmce_list = ET + FL + entity_ids = {} + cur = 0 -def e_lambda_1_2d(i, dx, dy, x_mid, y_mid): - EL = tuple([(0, y_mid**j*a) for a in dx for j in range(i)] + [(x_mid**j*a, 0) for a in dy for j in range(i)]) + for top_dim, entities in flat_topology.items(): + entity_ids[top_dim] = {} + for entity in entities: + entity_ids[top_dim][entity] = [] + + for j in sorted(flat_topology[1]): + entity_ids[1][j] = list(range(cur, cur + degree + 1)) + cur = cur + degree + 1 + + entity_ids[2][0] = list(range(cur, cur + len(FL))) + cur += len(FL) + + assert len(bdmce_list) == cur + + +def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): + EL = tuple([(0, y_mid**j*dx[0]) for j in range(deg)] + + [(y_mid**(r-1)*dy[0]*dy[1], (r+1)*y_mid**r*dx[0])] + + [(0, y_mid**j*dx[1]) for j in range(deg)] + + [(y_mid**(r-1)*dy[0]*dy[1], (r+1)*y_mid**r*dx[1])] + + [(x_mid**j*dy[0], 0) for j in range(deg)] + + [((r+1)*x_mid**r*dy[0], x_mid**(r-1)*dx[0]*dx[1])] + + [(x_mid**j*dy[1], 0) for j in range(deg)] + + [((r+1)*x_mid**r*dy[1], x_mid**(r-1)*dx[0]*dx[1])]) return EL -def e_tilda_lambda_1_2d(r, dx, dy, x_mid, y_mid): - ETL = tuple([(y_mid**(r-1)*dy[0]*dy[1], (r+1)*y_mid**r*a) for a in dx] + - [((r+1)*x_mid**r*a, x_mid**(r-1)*dx[0]*dx[1]) for a in dy]) - return ETL +def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): + FL = tuple([(x_mid**j*y_mid**(deg-2-j)*dy[0]*dy[1], 0) for j in range(2, deg+1)] + + [(0, x_mid**j*y_mid**(deg-2-j)*dx[0]*dx[1]) for j in range(2, deg+1)]) + + return FL From 6658adac18602100bdbd658bf8aa09983a4205fb Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Wed, 26 Jun 2019 01:00:11 +0100 Subject: [PATCH 03/23] first draft bdmc --- FIAT/brezzi_douglas_marini_cube.py | 62 ++++++++++++++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index ce6cd9d85..fbb75ffce 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -69,6 +69,68 @@ def __init__(self, ref_el, degree): cur += len(FL) assert len(bdmce_list) == cur + formdegree = 1 + + super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree) + + self.basis = {(0, 0): Array(bdmce_list)} + + topology = ref_el.get_topology() + unflattening_map = compute_unflattening_map(topology) + unflattened_entity_ids = {} + unflattened_entity_closure_ids = {} + + entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) + + for dim, entities in sorted(topology.items()): + unflattened_entity_ids[dim] = {} + unflattened_entity_closure_ids[dim] = {} + for dim, entities in sorted(flat_topology.items()): + for entity in entities: + unflat_dim, unflat_entity = unflattening_map[(dim, entity)] + unflattened_entity_ids[unflat_dim][unflat_entity] = entity_ids[dim][entity] + unflattened_entity_closure_ids[unflat_dim][unflat_entity] = entity_closure_ids[dim][entity] + self.entity_ids = unflattened_entity_ids + self.entity_closure_ids = unflattened_entity_closure_ids + self._degree = degree + self.flat_el = flat_el + + def degree(self): + return self._degree + 1 + + def get_nodal_basis(self): + raise NotImplementedError("get_nodal_basis not implemented for serendipity") + + def get_dual_set(self): + raise NotImplementedError("get_dual_set is not implemented for serendipity") + + def get_coeffs(self): + raise NotImplementedError("get_coeffs not implemented for serendipity") + + def tabulate(self, order, points, entity=None): + raise NotImplementedError + + def entity_dofs(self): + """Return the map of topological entities to degrees of + freedom for the finite element.""" + return self.entity_ids + + def entity_closure_dofs(self): + """Return the map of topological entities to degrees of + freedom on the closure of those entities for the finite element.""" + return self.entity_closure_ids + + def value_shape(self): + return (2,) + + def dmats(self): + raise NotImplementedError + + def get_num_members(self, arg): + raise NotImplementedError + + def space_dimension(self): + return len(self.basis[(0,)*self.flat_el.get_spatial_dimension()]) def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): From a7c5edeb069579ec46e73ecc9d4d9d704f616bcb Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Wed, 26 Jun 2019 01:50:48 +0100 Subject: [PATCH 04/23] minor fixes --- FIAT/brezzi_douglas_marini_cube.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index fbb75ffce..3de8573c3 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -24,6 +24,7 @@ from FIAT.dual_set import make_entity_closure_ids from FIAT.polynomial_set import mis from FIAT.serendipity import tr +from FIAT.reference_element import compute_unflattening_map, flatten_reference_cube x, y, z = symbols('x y z') variables = (x, y, z) @@ -51,7 +52,7 @@ def __init__(self, ref_el, degree): EL = e_lambda_1_2d(degree, dx, dy, x_mid, y_mid) FL = f_lambda_1_2d(degree, dx, dy, x_mid, y_mid) - bdmce_list = ET + FL + bdmce_list = EL + FL entity_ids = {} cur = 0 @@ -99,13 +100,13 @@ def degree(self): return self._degree + 1 def get_nodal_basis(self): - raise NotImplementedError("get_nodal_basis not implemented for serendipity") + raise NotImplementedError("get_nodal_basis not implemented for bdmce") def get_dual_set(self): - raise NotImplementedError("get_dual_set is not implemented for serendipity") + raise NotImplementedError("get_dual_set is not implemented for bdmce") def get_coeffs(self): - raise NotImplementedError("get_coeffs not implemented for serendipity") + raise NotImplementedError("get_coeffs not implemented for bdmce") def tabulate(self, order, points, entity=None): raise NotImplementedError @@ -135,13 +136,13 @@ def space_dimension(self): def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): EL = tuple([(0, y_mid**j*dx[0]) for j in range(deg)] + - [(y_mid**(r-1)*dy[0]*dy[1], (r+1)*y_mid**r*dx[0])] + + [(y_mid**(deg-1)*dy[0]*dy[1], (deg+1)*y_mid**deg*dx[0])] + [(0, y_mid**j*dx[1]) for j in range(deg)] + - [(y_mid**(r-1)*dy[0]*dy[1], (r+1)*y_mid**r*dx[1])] + + [(y_mid**(deg-1)*dy[0]*dy[1], (deg+1)*y_mid**deg*dx[1])] + [(x_mid**j*dy[0], 0) for j in range(deg)] + - [((r+1)*x_mid**r*dy[0], x_mid**(r-1)*dx[0]*dx[1])] + + [((deg+1)*x_mid**deg*dy[0], x_mid**(deg-1)*dx[0]*dx[1])] + [(x_mid**j*dy[1], 0) for j in range(deg)] + - [((r+1)*x_mid**r*dy[1], x_mid**(r-1)*dx[0]*dx[1])]) + [((deg+1)*x_mid**deg*dy[1], x_mid**(deg-1)*dx[0]*dx[1])]) return EL From aff450f27ff61c4affa1fa095e8f5929acd220f8 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Wed, 26 Jun 2019 13:01:23 +0100 Subject: [PATCH 05/23] change name poly to polynomials --- FIAT/serendipity.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/FIAT/serendipity.py b/FIAT/serendipity.py index 18c0cc7fa..c5a33ab8d 100644 --- a/FIAT/serendipity.py +++ b/FIAT/serendipity.py @@ -22,7 +22,7 @@ from FIAT.finite_element import FiniteElement from FIAT.lagrange import Lagrange from FIAT.dual_set import make_entity_closure_ids -from FIAT.polynomial_set import mis +from FIAT.polynomialsnomial_set import mis from FIAT.reference_element import (compute_unflattening_map, flatten_reference_cube) @@ -165,14 +165,14 @@ def tabulate(self, order, points, entity=None): alphas = mis(dim, o) for alpha in alphas: try: - poly = self.basis[alpha] + polynomials = self.basis[alpha] except KeyError: - poly = diff(self.basis[(0,)*dim], *zip(variables, alpha)) - self.basis[alpha] = poly - T = np.zeros((len(poly), len(points))) + polynomials = diff(self.basis[(0,)*dim], *zip(variables, alpha)) + self.basis[alpha] = polynomials + T = np.zeros((len(polynomials), len(points))) for i in range(len(points)): subs = {v: points[i][k] for k, v in enumerate(variables[:dim])} - for j, f in enumerate(poly): + for j, f in enumerate(polynomials): T[j, i] = f.evalf(subs=subs) phivals[alpha] = T From 276512c430f75823d119907fd9e126b603e2199d Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Wed, 26 Jun 2019 14:11:31 +0100 Subject: [PATCH 06/23] minor updates --- FIAT/brezzi_douglas_marini_cube.py | 34 ++++++++++++++++++++++++++---- FIAT/reference_element.py | 23 ++++++++++++++------ FIAT/serendipity.py | 2 +- 3 files changed, 48 insertions(+), 11 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 3de8573c3..c1d60213f 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -72,6 +72,8 @@ def __init__(self, ref_el, degree): assert len(bdmce_list) == cur formdegree = 1 + entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) + super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree) self.basis = {(0, 0): Array(bdmce_list)} @@ -81,8 +83,6 @@ def __init__(self, ref_el, degree): unflattened_entity_ids = {} unflattened_entity_closure_ids = {} - entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) - for dim, entities in sorted(topology.items()): unflattened_entity_ids[dim] = {} unflattened_entity_closure_ids[dim] = {} @@ -109,7 +109,33 @@ def get_coeffs(self): raise NotImplementedError("get_coeffs not implemented for bdmce") def tabulate(self, order, points, entity=None): - raise NotImplementedError + + if entity is None: + entity = (self.ref_el.get_spatial_dimension(), 0) + + entity_dim, entity_id = entity + transform = self.ref_el.get_entity_transform(entity_dim, entity_id) + points = list(map(transform, points)) + + phivals = {} + + for o in range(order+1): + alphas = mis(2, o) + for alpha in alphas: + try: + polynomials = self.basis[alpha] + except KeyError: + polynomials = diff(self.basis[(0,)*dim], *zip(variables, alpha)) + self.basis[alpha] = polynomials + T = np.zeros((len(polynomials), 2, len(points))) + for i in range(len(points)): + subs = {v: points[i][k] for k, v in enumerate(variables[:2])} + for k in range(2): + for j, f in enumerate(polynomials[:, k]): + T[j, :, i] = f.evalf(subs=subs) + phivals[alpha] = T + + return phivals def entity_dofs(self): """Return the map of topological entities to degrees of @@ -131,7 +157,7 @@ def get_num_members(self, arg): raise NotImplementedError def space_dimension(self): - return len(self.basis[(0,)*self.flat_el.get_spatial_dimension()]) + return 2 def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index e2c38f191..259546e3d 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1025,16 +1025,27 @@ def tuple_sum(tree): return tree +def is_hypercube(cell): + if isinstance(cell, (DefaultLine, UFCInterval, UFCQuadrilateral, UFCHexahedron)): + return True + elif isinstance(cell, TensorProductCell): + return reduce(lambda a, b: a and b, [is_hypercube(c) for c in cell.cells]) + else: + return False + def flatten_reference_cube(ref_el): - """Returns a flattened cube corresponding to a given tensor product cell.""" + """ ... """ + flattened_cube = {2: UFCQuadrilateral(), 3: UFCHexahedron()} if numpy.sum(ref_el.get_dimension()) <= 1: + # Just return point/interval cell arguments return ref_el - elif numpy.sum(ref_el.get_dimension()) == 2 and len(ref_el.vertices) == 4: - return UFCQuadrilateral() - elif numpy.sum(ref_el.get_dimension()) == 3 and len(ref_el.vertices) == 8: - return UFCHexahedron() else: - raise TypeError('Invalid cell type') + # Handle cases where cell is a quad/cube constructed from a tensor product or + # an already flattened element + if is_hypercube(ref_el): + return flattened_cube[numpy.sum(ref_el.get_dimension())] + else: + raise TypeError('Invalid cell type') def flatten_entities(topology_dict): diff --git a/FIAT/serendipity.py b/FIAT/serendipity.py index c5a33ab8d..ab0cac1a2 100644 --- a/FIAT/serendipity.py +++ b/FIAT/serendipity.py @@ -22,7 +22,7 @@ from FIAT.finite_element import FiniteElement from FIAT.lagrange import Lagrange from FIAT.dual_set import make_entity_closure_ids -from FIAT.polynomialsnomial_set import mis +from FIAT.polynomial_set import mis from FIAT.reference_element import (compute_unflattening_map, flatten_reference_cube) From de2976c8d2288d5be8b13e71619e564437558cf1 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Wed, 26 Jun 2019 14:54:35 +0100 Subject: [PATCH 07/23] flake8 complaint --- FIAT/reference_element.py | 1 + 1 file changed, 1 insertion(+) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 259546e3d..76ef111eb 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1033,6 +1033,7 @@ def is_hypercube(cell): else: return False + def flatten_reference_cube(ref_el): """ ... """ flattened_cube = {2: UFCQuadrilateral(), 3: UFCHexahedron()} From 7ebba5a81f08f97c841612d6ee7ae903315f3ab2 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Thu, 27 Jun 2019 21:58:49 +0100 Subject: [PATCH 08/23] minor edits --- FIAT/__init__.py | 2 ++ FIAT/brezzi_douglas_marini_cube.py | 31 +++++++++++++++++------------- 2 files changed, 20 insertions(+), 13 deletions(-) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index f5c1677ce..59f72ce7f 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -17,6 +17,7 @@ from FIAT.discontinuous_raviart_thomas import DiscontinuousRaviartThomas from FIAT.serendipity import Serendipity from FIAT.discontinuous_pc import DPC +from FIAT.brezzi_douglas_marini_cube import BrezziDouglasMariniCube from FIAT.hermite import CubicHermite from FIAT.lagrange import Lagrange from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre @@ -58,6 +59,7 @@ "Crouzeix-Raviart": CrouzeixRaviart, "Discontinuous Lagrange": DiscontinuousLagrange, "S": Serendipity, + "Brezzi-Douglas-Marini-Cube": BrezziDouglasMariniCube, "DPC": DPC, "Discontinuous Taylor": DiscontinuousTaylor, "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas, diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index c1d60213f..552c786ab 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -31,7 +31,7 @@ leg = legendre -class BrezziDouglasMariniCubeEdge(FiniteElement): +class BrezziDouglasMariniCube(FiniteElement): def __init__(self, ref_el, degree): if degree < 1: raise Exception("BDMce_k elements only valid for k >= 1") @@ -51,7 +51,10 @@ def __init__(self, ref_el, degree): y_mid = 2*y-(verts[-1][1] + verts[0][1]) EL = e_lambda_1_2d(degree, dx, dy, x_mid, y_mid) - FL = f_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + if degree >= 2: + FL = f_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + else: + FL = () bdmce_list = EL + FL entity_ids = {} @@ -74,7 +77,8 @@ def __init__(self, ref_el, degree): entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) - super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree) + super(BrezziDouglasMariniCube, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, + mapping="contravariant piola") self.basis = {(0, 0): Array(bdmce_list)} @@ -97,7 +101,7 @@ def __init__(self, ref_el, degree): self.flat_el = flat_el def degree(self): - return self._degree + 1 + return self._degree def get_nodal_basis(self): raise NotImplementedError("get_nodal_basis not implemented for bdmce") @@ -111,7 +115,7 @@ def get_coeffs(self): def tabulate(self, order, points, entity=None): if entity is None: - entity = (self.ref_el.get_spatial_dimension(), 0) + entity = (self.ref_el.get_dimension(), 0) entity_dim, entity_id = entity transform = self.ref_el.get_entity_transform(entity_dim, entity_id) @@ -125,14 +129,15 @@ def tabulate(self, order, points, entity=None): try: polynomials = self.basis[alpha] except KeyError: - polynomials = diff(self.basis[(0,)*dim], *zip(variables, alpha)) + polynomials = diff(self.basis[(0, 0)], *zip(variables, alpha)) self.basis[alpha] = polynomials - T = np.zeros((len(polynomials), 2, len(points))) + T = np.zeros((len(polynomials[:, 0]), 2, len(points))) for i in range(len(points)): subs = {v: points[i][k] for k, v in enumerate(variables[:2])} - for k in range(2): - for j, f in enumerate(polynomials[:, k]): - T[j, :, i] = f.evalf(subs=subs) + for j, f in enumerate(polynomials[:, 0]): + T[j, 0, i] = f.evalf(subs=subs) + for j, f in enumerate(polynomials[:, 1]): + T[j, 1, i] = f.evalf(subs=subs) phivals[alpha] = T return phivals @@ -157,7 +162,7 @@ def get_num_members(self, arg): raise NotImplementedError def space_dimension(self): - return 2 + return int(len(self.basis[(0, 0)])/2) def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): @@ -174,7 +179,7 @@ def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): - FL = tuple([(x_mid**j*y_mid**(deg-2-j)*dy[0]*dy[1], 0) for j in range(2, deg+1)] + - [(0, x_mid**j*y_mid**(deg-2-j)*dx[0]*dx[1]) for j in range(2, deg+1)]) + FL = tuple([(x_mid**j*y_mid**(deg-2-j)*dy[0]*dy[1], 0) for j in range(deg-1)] + + [(0, x_mid**j*y_mid**(deg-2-j)*dx[0]*dx[1]) for j in range(deg-1)]) return FL From 13b6abba87016285bc5a9ed1d90c88acf05390af Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 28 Jun 2019 00:01:34 +0100 Subject: [PATCH 09/23] minor changes --- FIAT/__init__.py | 1 + FIAT/brezzi_douglas_marini_cube.py | 36 ++++++++++++++++++++++++++++++ 2 files changed, 37 insertions(+) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 59f72ce7f..e3d08cc7b 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -16,6 +16,7 @@ from FIAT.discontinuous_taylor import DiscontinuousTaylor from FIAT.discontinuous_raviart_thomas import DiscontinuousRaviartThomas from FIAT.serendipity import Serendipity +from FIAT.brezzi_douglas_marini_cube import BrezziDouglasMariniCubeEdge, BrezziDouglasMariniCubeFace from FIAT.discontinuous_pc import DPC from FIAT.brezzi_douglas_marini_cube import BrezziDouglasMariniCube from FIAT.hermite import CubicHermite diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 552c786ab..661ae17bc 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -183,3 +183,39 @@ def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): [(0, x_mid**j*y_mid**(deg-2-j)*dx[0]*dx[1]) for j in range(deg-1)]) return FL + + +class BrezziDouglasMariniCubeEdge(BrezziDouglasMariniCube): + def __init__(self, ref_el, degree): + super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, degree=degree) + + +class BrezziDouglasMariniCubeFace(BrezziDouglasMariniCube): + def __init__(self, ref_el, degree): + if degree < 1: + raise Exception("BDMcf_k elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + if dim != 2: + raise Exception("BDMcf_k elements only valid for dimension 2") + + flat_topology = flat_el.get_topology() + + verts = flat_el.get_vertices() + + dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) + dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) + x_mid = 2*x-(verts[-1][0] + verts[0][0]) + y_mid = 2*y-(verts[-1][1] + verts[0][1]) + + EL = e_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + if degree >= 2: + FL = f_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + else: + FL = () + bdmcf_list = EL + FL + bdmcf_list = [[a[1], -a[0]] for a in bdmcf_list] + self.basis = {(0, 0): Array(bdmcf_list)} + + super(BrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree) From ab2848c858ef12b6bfddeb3a92a23fdb697d53ad Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 28 Jun 2019 19:57:25 +0100 Subject: [PATCH 10/23] minor changes --- FIAT/brezzi_douglas_marini_cube.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 661ae17bc..d71254e48 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -166,14 +166,14 @@ def space_dimension(self): def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): - EL = tuple([(0, y_mid**j*dx[0]) for j in range(deg)] + - [(y_mid**(deg-1)*dy[0]*dy[1], (deg+1)*y_mid**deg*dx[0])] + - [(0, y_mid**j*dx[1]) for j in range(deg)] + - [(y_mid**(deg-1)*dy[0]*dy[1], (deg+1)*y_mid**deg*dx[1])] + - [(x_mid**j*dy[0], 0) for j in range(deg)] + - [((deg+1)*x_mid**deg*dy[0], x_mid**(deg-1)*dx[0]*dx[1])] + - [(x_mid**j*dy[1], 0) for j in range(deg)] + - [((deg+1)*x_mid**deg*dy[1], x_mid**(deg-1)*dx[0]*dx[1])]) + EL = tuple([(0, -y_mid**j*dx[0]) for j in range(deg)] + + [(-y_mid**(deg-1)*dy[0]*dy[1], -(deg+1)*y_mid**deg*dx[0])] + + [(0, -y_mid**j*dx[1]) for j in range(deg)] + + [(-y_mid**(deg-1)*dy[0]*dy[1], -(deg+1)*y_mid**deg*dx[1])] + + [(-x_mid**j*dy[0], 0) for j in range(deg)] + + [(-(deg+1)*x_mid**deg*dy[0], -x_mid**(deg-1)*dx[0]*dx[1])] + + [(-x_mid**j*dy[1], 0) for j in range(deg)] + + [(-(deg+1)*x_mid**deg*dy[1], -x_mid**(deg-1)*dx[0]*dx[1])]) return EL @@ -217,5 +217,5 @@ def __init__(self, ref_el, degree): bdmcf_list = EL + FL bdmcf_list = [[a[1], -a[0]] for a in bdmcf_list] self.basis = {(0, 0): Array(bdmcf_list)} - + super(BrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree) From 78d9395291862d40c50b8e0223887bc0eb087e93 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 12 Jul 2019 11:31:20 +0100 Subject: [PATCH 11/23] new docstring for flatten_reference_cube --- FIAT/reference_element.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 76ef111eb..232cbf767 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1035,7 +1035,7 @@ def is_hypercube(cell): def flatten_reference_cube(ref_el): - """ ... """ + """This function flattens a Tensor Product hypercube to the corresponding UFC hypercube""" flattened_cube = {2: UFCQuadrilateral(), 3: UFCHexahedron()} if numpy.sum(ref_el.get_dimension()) <= 1: # Just return point/interval cell arguments From cd03c26effadad2b63ce6db6a2ad8f7a51664f4e Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 12 Jul 2019 14:44:03 +0100 Subject: [PATCH 12/23] fix docstring --- FIAT/serendipity.py | 1 - 1 file changed, 1 deletion(-) diff --git a/FIAT/serendipity.py b/FIAT/serendipity.py index ab0cac1a2..eca81442e 100644 --- a/FIAT/serendipity.py +++ b/FIAT/serendipity.py @@ -56,7 +56,6 @@ def __init__(self, ref_el, degree): dim = flat_el.get_spatial_dimension() flat_topology = flat_el.get_topology() - x, y, z = symbols('x y z') verts = flat_el.get_vertices() dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) From 87841f9d74b8e22246bf672548c13fbd625d2f6d Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Sat, 27 Jul 2019 20:52:58 +0100 Subject: [PATCH 13/23] changed basis function order --- FIAT/brezzi_douglas_marini_cube.py | 71 +++++++++++++++++------------- 1 file changed, 41 insertions(+), 30 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index d71254e48..683b43721 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -31,6 +31,10 @@ leg = legendre +def triangular_number(n): + return int((n+1)*n/2) + + class BrezziDouglasMariniCube(FiniteElement): def __init__(self, ref_el, degree): if degree < 1: @@ -43,20 +47,6 @@ def __init__(self, ref_el, degree): flat_topology = flat_el.get_topology() - verts = flat_el.get_vertices() - - dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) - dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) - x_mid = 2*x-(verts[-1][0] + verts[0][0]) - y_mid = 2*y-(verts[-1][1] + verts[0][1]) - - EL = e_lambda_1_2d(degree, dx, dy, x_mid, y_mid) - if degree >= 2: - FL = f_lambda_1_2d(degree, dx, dy, x_mid, y_mid) - else: - FL = () - bdmce_list = EL + FL - entity_ids = {} cur = 0 @@ -69,10 +59,9 @@ def __init__(self, ref_el, degree): entity_ids[1][j] = list(range(cur, cur + degree + 1)) cur = cur + degree + 1 - entity_ids[2][0] = list(range(cur, cur + len(FL))) - cur += len(FL) + entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 1))) + cur += 2*triangular_number(degree - 1) - assert len(bdmce_list) == cur formdegree = 1 entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) @@ -80,8 +69,6 @@ def __init__(self, ref_el, degree): super(BrezziDouglasMariniCube, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, mapping="contravariant piola") - self.basis = {(0, 0): Array(bdmce_list)} - topology = ref_el.get_topology() unflattening_map = compute_unflattening_map(topology) unflattened_entity_ids = {} @@ -166,27 +153,51 @@ def space_dimension(self): def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): - EL = tuple([(0, -y_mid**j*dx[0]) for j in range(deg)] + - [(-y_mid**(deg-1)*dy[0]*dy[1], -(deg+1)*y_mid**deg*dx[0])] + - [(0, -y_mid**j*dx[1]) for j in range(deg)] + - [(-y_mid**(deg-1)*dy[0]*dy[1], -(deg+1)*y_mid**deg*dx[1])] + - [(-x_mid**j*dy[0], 0) for j in range(deg)] + - [(-(deg+1)*x_mid**deg*dy[0], -x_mid**(deg-1)*dx[0]*dx[1])] + - [(-x_mid**j*dy[1], 0) for j in range(deg)] + - [(-(deg+1)*x_mid**deg*dy[1], -x_mid**(deg-1)*dx[0]*dx[1])]) + EL = tuple([(0, -leg(j, y_mid)*dx[0]) for j in range(deg)] + + [(-leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[0])] + + [(0, -leg(j, y_mid)*dx[1]) for j in range(deg)] + + [(leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[1])] + + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + + [(-leg(deg, x_mid)*dy[0], -leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))] + + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)] + + [(-leg(deg, x_mid)*dy[1], leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))]) return EL def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): - FL = tuple([(x_mid**j*y_mid**(deg-2-j)*dy[0]*dy[1], 0) for j in range(deg-1)] + - [(0, x_mid**j*y_mid**(deg-2-j)*dx[0]*dx[1]) for j in range(deg-1)]) + FL = tuple([(x_mid**j*y_mid**(k-2-j)*dy[0]*dy[1], 0) for k in range(2, deg+1) for j in range(k-1)] + + [(0, x_mid**j*y_mid**(k-2-j)*dx[0]*dx[1]) for k in range(2, deg+1) for j in range(k-1)]) return FL class BrezziDouglasMariniCubeEdge(BrezziDouglasMariniCube): def __init__(self, ref_el, degree): + if degree < 1: + raise Exception("BDMcf_k elements only valid for k >= 1") + + flat_el = flatten_reference_cube(ref_el) + dim = flat_el.get_spatial_dimension() + if dim != 2: + raise Exception("BDMcf_k elements only valid for dimension 2") + + flat_topology = flat_el.get_topology() + + verts = flat_el.get_vertices() + + dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) + dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) + x_mid = 2*x-(verts[-1][0] + verts[0][0]) + y_mid = 2*y-(verts[-1][1] + verts[0][1]) + + EL = e_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + if degree >= 2: + FL = f_lambda_1_2d(degree, dx, dy, x_mid, y_mid) + else: + FL = () + bdmce_list = EL + FL + self.basis = {(0, 0): Array(bdmce_list)} super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, degree=degree) @@ -215,7 +226,7 @@ def __init__(self, ref_el, degree): else: FL = () bdmcf_list = EL + FL - bdmcf_list = [[a[1], -a[0]] for a in bdmcf_list] + bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] self.basis = {(0, 0): Array(bdmcf_list)} super(BrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree) From 2c42a0ca8a344bac14ffb22f4a5ff57c26c72636 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Thu, 8 Aug 2019 01:31:36 +0100 Subject: [PATCH 14/23] working for degree 1 and 2 --- FIAT/brezzi_douglas_marini_cube.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 683b43721..96cf2afdf 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -154,20 +154,20 @@ def space_dimension(self): def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): EL = tuple([(0, -leg(j, y_mid)*dx[0]) for j in range(deg)] + - [(-leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[0])] + + [(-leg(deg-1, y_mid)*dy[0]*dy[1]/deg, -leg(deg, y_mid)*dx[0])] + [(0, -leg(j, y_mid)*dx[1]) for j in range(deg)] + - [(leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[1])] + + [(leg(deg-1, y_mid)*dy[0]*dy[1]/deg, -leg(deg, y_mid)*dx[1])] + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + - [(-leg(deg, x_mid)*dy[0], -leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))] + + [(-leg(deg, x_mid)*dy[0], -leg(deg-1, x_mid)*dx[0]*dx[1]/deg)] + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)] + - [(-leg(deg, x_mid)*dy[1], leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))]) + [(-leg(deg, x_mid)*dy[1], leg(deg-1, x_mid)*dx[0]*dx[1]/deg)]) return EL def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): - FL = tuple([(x_mid**j*y_mid**(k-2-j)*dy[0]*dy[1], 0) for k in range(2, deg+1) for j in range(k-1)] + - [(0, x_mid**j*y_mid**(k-2-j)*dx[0]*dx[1]) for k in range(2, deg+1) for j in range(k-1)]) + FL = tuple([(leg(j, x_mid)*leg(k-2-j, y_mid)*dy[0]*dy[1], 0) for k in range(2, deg+1) for j in range(k-1)] + + [(0, -leg(k-2-j, x_mid)*leg(j, y_mid)*dx[0]*dx[1]) for k in range(2, deg+1) for j in range(k-1)]) return FL From c84d83d64d7f62b6c178fa37d58cf2e12677fcbf Mon Sep 17 00:00:00 2001 From: David Ham Date: Thu, 8 Aug 2019 17:23:26 +0100 Subject: [PATCH 15/23] Flake8 --- FIAT/__init__.py | 4 ++-- FIAT/brezzi_douglas_marini_cube.py | 6 ------ 2 files changed, 2 insertions(+), 8 deletions(-) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index e3d08cc7b..327669b98 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -18,7 +18,6 @@ from FIAT.serendipity import Serendipity from FIAT.brezzi_douglas_marini_cube import BrezziDouglasMariniCubeEdge, BrezziDouglasMariniCubeFace from FIAT.discontinuous_pc import DPC -from FIAT.brezzi_douglas_marini_cube import BrezziDouglasMariniCube from FIAT.hermite import CubicHermite from FIAT.lagrange import Lagrange from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre @@ -60,7 +59,8 @@ "Crouzeix-Raviart": CrouzeixRaviart, "Discontinuous Lagrange": DiscontinuousLagrange, "S": Serendipity, - "Brezzi-Douglas-Marini-Cube": BrezziDouglasMariniCube, + "Brezzi-Douglas-Marini Cube Face": BrezziDouglasMariniCubeFace, + "Brezzi-Douglas-Marini Cube Edge": BrezziDouglasMariniCubeEdge, "DPC": DPC, "Discontinuous Taylor": DiscontinuousTaylor, "Discontinuous Raviart-Thomas": DiscontinuousRaviartThomas, diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 96cf2afdf..8d66d5716 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -20,10 +20,8 @@ from sympy import symbols, legendre, Array, diff import numpy as np from FIAT.finite_element import FiniteElement -from FIAT.lagrange import Lagrange from FIAT.dual_set import make_entity_closure_ids from FIAT.polynomial_set import mis -from FIAT.serendipity import tr from FIAT.reference_element import compute_unflattening_map, flatten_reference_cube x, y, z = symbols('x y z') @@ -182,8 +180,6 @@ def __init__(self, ref_el, degree): if dim != 2: raise Exception("BDMcf_k elements only valid for dimension 2") - flat_topology = flat_el.get_topology() - verts = flat_el.get_vertices() dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) @@ -211,8 +207,6 @@ def __init__(self, ref_el, degree): if dim != 2: raise Exception("BDMcf_k elements only valid for dimension 2") - flat_topology = flat_el.get_topology() - verts = flat_el.get_vertices() dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) From d33bac146175d63fda2d5103e6f0982f8fa1275f Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 9 Aug 2019 03:03:50 +0100 Subject: [PATCH 16/23] working bdmcf --- FIAT/brezzi_douglas_marini_cube.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 8d66d5716..a14432d0a 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -164,10 +164,13 @@ def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): - FL = tuple([(leg(j, x_mid)*leg(k-2-j, y_mid)*dy[0]*dy[1], 0) for k in range(2, deg+1) for j in range(k-1)] + - [(0, -leg(k-2-j, x_mid)*leg(j, y_mid)*dx[0]*dx[1]) for k in range(2, deg+1) for j in range(k-1)]) + FL = [] + for k in range(2, deg+1): + for j in range(k-1): + FL += [(0, leg(j, x_mid)*leg(k-2-j, y_mid)*dx[0]*dx[1])] + FL += [(leg(k-2-j, x_mid)*leg(j, y_mid)*dy[0]*dy[1], 0)] - return FL + return tuple(FL) class BrezziDouglasMariniCubeEdge(BrezziDouglasMariniCube): From 48b59d21e32f25b92053f36f01df6093caf2640f Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Fri, 9 Aug 2019 14:52:46 +0100 Subject: [PATCH 17/23] Andrew's basis functions --- FIAT/brezzi_douglas_marini_cube.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index a14432d0a..942a9c017 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -152,13 +152,13 @@ def space_dimension(self): def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): EL = tuple([(0, -leg(j, y_mid)*dx[0]) for j in range(deg)] + - [(-leg(deg-1, y_mid)*dy[0]*dy[1]/deg, -leg(deg, y_mid)*dx[0])] + + [(-leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[0])] + [(0, -leg(j, y_mid)*dx[1]) for j in range(deg)] + - [(leg(deg-1, y_mid)*dy[0]*dy[1]/deg, -leg(deg, y_mid)*dx[1])] + + [(leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[1])] + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + - [(-leg(deg, x_mid)*dy[0], -leg(deg-1, x_mid)*dx[0]*dx[1]/deg)] + + [(-leg(deg, x_mid)*dy[0], -leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))] + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)] + - [(-leg(deg, x_mid)*dy[1], leg(deg-1, x_mid)*dx[0]*dx[1]/deg)]) + [(-leg(deg, x_mid)*dy[1], leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))]) return EL From a4f88e4e1dff28adf2fbbfc51ad4f279e1e0d0c0 Mon Sep 17 00:00:00 2001 From: cyruscycheng21 Date: Sat, 24 Aug 2019 18:17:22 +0100 Subject: [PATCH 18/23] add mappings --- FIAT/brezzi_douglas_marini_cube.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 942a9c017..44856aac5 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -34,7 +34,7 @@ def triangular_number(n): class BrezziDouglasMariniCube(FiniteElement): - def __init__(self, ref_el, degree): + def __init__(self, ref_el, degree, mapping): if degree < 1: raise Exception("BDMce_k elements only valid for k >= 1") @@ -65,7 +65,7 @@ def __init__(self, ref_el, degree): entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) super(BrezziDouglasMariniCube, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, - mapping="contravariant piola") + mapping=mapping) topology = ref_el.get_topology() unflattening_map = compute_unflattening_map(topology) @@ -197,7 +197,7 @@ def __init__(self, ref_el, degree): FL = () bdmce_list = EL + FL self.basis = {(0, 0): Array(bdmce_list)} - super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, degree=degree) + super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") class BrezziDouglasMariniCubeFace(BrezziDouglasMariniCube): @@ -226,4 +226,4 @@ def __init__(self, ref_el, degree): bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] self.basis = {(0, 0): Array(bdmcf_list)} - super(BrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree) + super(BrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") From 758ef376fae4cc161b0e367d3ccf56e2ea053fd8 Mon Sep 17 00:00:00 2001 From: Tom Bendall Date: Tue, 26 Jan 2021 19:17:03 +0000 Subject: [PATCH 19/23] fix BDMCE/F elements by correcting coefficients for basis --- FIAT/brezzi_douglas_marini_cube.py | 24 +++++++++++++++++++----- 1 file changed, 19 insertions(+), 5 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 44856aac5..38285fcdb 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -17,7 +17,7 @@ # # Modified by David A. Ham (david.ham@imperial.ac.uk), 2019 -from sympy import symbols, legendre, Array, diff +from sympy import symbols, legendre, Array, diff, binomial import numpy as np from FIAT.finite_element import FiniteElement from FIAT.dual_set import make_entity_closure_ids @@ -151,14 +151,28 @@ def space_dimension(self): def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): + # Elements introduced by Brezzi, Douglas, Marini (1985) + # "Two families of mixed finite elements for Second Order Elliptic Problems" + + # See e.g. Brezzi, Douglas, Fortin, Marini (1987) + # "Efficient rectangular mixed finite elements in two and three space variables" + # For rectangle K: + # BDM_j(K) = [P_k(K)^2 + Span(curl(xy^{j+1}, x^{j+1}y))] x P_{j-1}(K) + # + # For all basis functions, the div/curl should be of degree (deg-1) + # Thus the div/curl component needs weighting by a coefficient + # This comes from looking at the coefficient of the highest order polynomial + # For leg(deg, 2x-1), this is binomial(2*deg, deg) + # The then takes this factor from both components of the vector + coeff = binomial(2*deg, deg) / ((deg+1)*binomial(2*deg-2, deg-1)) EL = tuple([(0, -leg(j, y_mid)*dx[0]) for j in range(deg)] + - [(-leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[0])] + + [(coeff*-leg(deg-1, y_mid)*dy[0]*dy[1], -leg(deg, y_mid)*dx[0])] + [(0, -leg(j, y_mid)*dx[1]) for j in range(deg)] + - [(leg(deg-1, y_mid)*dy[0]*dy[1]/(deg+1), -leg(deg, y_mid)*dx[1])] + + [(coeff*leg(deg-1, y_mid)*dy[0]*dy[1], -leg(deg, y_mid)*dx[1])] + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + - [(-leg(deg, x_mid)*dy[0], -leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))] + + [(-leg(deg, x_mid)*dy[0], coeff*-leg(deg-1, x_mid)*dx[0]*dx[1])] + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)] + - [(-leg(deg, x_mid)*dy[1], leg(deg-1, x_mid)*dx[0]*dx[1]/(deg+1))]) + [(-leg(deg, x_mid)*dy[1], coeff*leg(deg-1, x_mid)*dx[0]*dx[1])]) return EL From d10132a2703285bf7fb357051e47d1a89d69422c Mon Sep 17 00:00:00 2001 From: Tom Bendall Date: Thu, 28 Jan 2021 11:19:11 +0000 Subject: [PATCH 20/23] Tidy up BDMCE/F elements --- FIAT/brezzi_douglas_marini_cube.py | 307 +++++++++++++++++++---------- 1 file changed, 206 insertions(+), 101 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 38285fcdb..efc6878c5 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -16,16 +16,17 @@ # along with FIAT. If not, see . # # Modified by David A. Ham (david.ham@imperial.ac.uk), 2019 +# Modified by Thomas Bendall (thomas.bendall@metoffice.gov.uk) 2021 -from sympy import symbols, legendre, Array, diff, binomial +from sympy import symbols, legendre, Array, diff, binomial, lambdify import numpy as np -from FIAT.finite_element import FiniteElement from FIAT.dual_set import make_entity_closure_ids +from FIAT.finite_element import FiniteElement from FIAT.polynomial_set import mis from FIAT.reference_element import compute_unflattening_map, flatten_reference_cube -x, y, z = symbols('x y z') -variables = (x, y, z) +x, y = symbols('x y') +variables = (x, y) leg = legendre @@ -34,19 +35,31 @@ def triangular_number(n): class BrezziDouglasMariniCube(FiniteElement): + """ + The Brezzi-Douglas-Marini element on quadrilateral cells. + + :arg ref_el: The reference element. + :arg k: The degree. + :arg mapping: A string giving the Piola mapping. + Either 'contravariant Piola' or 'covariant Piola'. + """ + def __init__(self, ref_el, degree, mapping): + + # Check that ref_el and degree are appropriate if degree < 1: - raise Exception("BDMce_k elements only valid for k >= 1") + raise Exception("BDMc_k elements only valid for k >= 1") flat_el = flatten_reference_cube(ref_el) dim = flat_el.get_spatial_dimension() if dim != 2: - raise Exception("BDMce_k elements only valid for dimension 2") + raise Exception("BDMc_k elements only valid for dimension 2") + # Collect the IDs of the reference element entities flat_topology = flat_el.get_topology() entity_ids = {} - cur = 0 + counter = 0 for top_dim, entities in flat_topology.items(): entity_ids[top_dim] = {} @@ -54,19 +67,20 @@ def __init__(self, ref_el, degree, mapping): entity_ids[top_dim][entity] = [] for j in sorted(flat_topology[1]): - entity_ids[1][j] = list(range(cur, cur + degree + 1)) - cur = cur + degree + 1 + entity_ids[1][j] = list(range(counter, counter + degree + 1)) + counter += degree + 1 - entity_ids[2][0] = list(range(cur, cur + 2*triangular_number(degree - 1))) - cur += 2*triangular_number(degree - 1) - - formdegree = 1 + entity_ids[2][0] = list(range(counter, counter + 2*triangular_number(degree - 1))) + counter += 2*triangular_number(degree - 1) entity_closure_ids = make_entity_closure_ids(flat_el, entity_ids) - super(BrezziDouglasMariniCube, self).__init__(ref_el=ref_el, dual=None, order=degree, formdegree=formdegree, + # Set up FiniteElement + super(BrezziDouglasMariniCube, self).__init__(ref_el=ref_el, dual=None, + order=degree, formdegree=1, mapping=mapping) + # Store unflattened entity ID dictionaries topology = ref_el.get_topology() unflattening_map = compute_unflattening_map(topology) unflattened_entity_ids = {} @@ -85,44 +99,63 @@ def __init__(self, ref_el, degree, mapping): self._degree = degree self.flat_el = flat_el + def degree(self): + """Return the degree of the polynomial space.""" return self._degree def get_nodal_basis(self): - raise NotImplementedError("get_nodal_basis not implemented for bdmce") + raise NotImplementedError("get_nodal_basis not implemented for BDMCE/F elements") def get_dual_set(self): - raise NotImplementedError("get_dual_set is not implemented for bdmce") + raise NotImplementedError("get_dual_set is not implemented for BDMCE/F elements") def get_coeffs(self): - raise NotImplementedError("get_coeffs not implemented for bdmce") + raise NotImplementedError("get_coeffs not implemented for BDMCE/F elements") def tabulate(self, order, points, entity=None): - + """Return tabulated values of derivatives up to a given order of + basis functions at given points. + + :arg order: The maximum order of derivative. + :arg points: An iterable of points. + :arg entity: Optional (dimension, entity number) pair + indicating which topological entity of the + reference element to tabulate on. If ``None``, + tabulated values are computed by geometrically + approximating which facet the points are on. + """ if entity is None: entity = (self.ref_el.get_dimension(), 0) entity_dim, entity_id = entity transform = self.ref_el.get_entity_transform(entity_dim, entity_id) - points = list(map(transform, points)) + points = np.asarray(points) + npoints, pointdim = points.shape - phivals = {} + # Turn analytic basis functions into python functions via lambdify + basis_callable = {(0, 0): numpy_lambdify(variables, self.basis[(0,0)], + modules="numpy")} + # Dictionary of values of functions + phivals = {} for o in range(order+1): alphas = mis(2, o) + # Collect basis and derivatives of basis for alpha in alphas: try: - polynomials = self.basis[alpha] + callable = basis_callable[alpha] except KeyError: - polynomials = diff(self.basis[(0, 0)], *zip(variables, alpha)) - self.basis[alpha] = polynomials - T = np.zeros((len(polynomials[:, 0]), 2, len(points))) - for i in range(len(points)): - subs = {v: points[i][k] for k, v in enumerate(variables[:2])} - for j, f in enumerate(polynomials[:, 0]): - T[j, 0, i] = f.evalf(subs=subs) - for j, f in enumerate(polynomials[:, 1]): - T[j, 1, i] = f.evalf(subs=subs) + diff_basis = diff(self.basis[(0,0)], *zip(variables, alpha)) + callable = numpy_lambdify(variables, diff_basis, modules="numpy") + self.basis[alpha] = diff_basis + + # tabulate by passing points through all lambdified functions + # resulting array has shape (len(self.basis), spatial_dim, npoints) + T = np.array([[[func_component(point) for point in points] + for func_component in func] + for func in callable]) + phivals[alpha] = T return phivals @@ -138,7 +171,8 @@ def entity_closure_dofs(self): return self.entity_closure_ids def value_shape(self): - return (2,) + """Return the value shape of the finite element functions.""" + return np.shape(self.basis[(0,0)][0]) def dmats(self): raise NotImplementedError @@ -147,97 +181,168 @@ def get_num_members(self, arg): raise NotImplementedError def space_dimension(self): + """Return the dimension of the finite element space.""" return int(len(self.basis[(0, 0)])/2) -def e_lambda_1_2d(deg, dx, dy, x_mid, y_mid): - # Elements introduced by Brezzi, Douglas, Marini (1985) - # "Two families of mixed finite elements for Second Order Elliptic Problems" - - # See e.g. Brezzi, Douglas, Fortin, Marini (1987) - # "Efficient rectangular mixed finite elements in two and three space variables" - # For rectangle K: - # BDM_j(K) = [P_k(K)^2 + Span(curl(xy^{j+1}, x^{j+1}y))] x P_{j-1}(K) - # - # For all basis functions, the div/curl should be of degree (deg-1) - # Thus the div/curl component needs weighting by a coefficient - # This comes from looking at the coefficient of the highest order polynomial - # For leg(deg, 2x-1), this is binomial(2*deg, deg) - # The then takes this factor from both components of the vector +def bdmce_edge_basis(deg, dx, dy, x_mid, y_mid): + """ Returns the basis functions associated with DoFs on the + edges of elements for the HCurl Brezz-Douglas-Marini element on + quadrilateral cells. + + These were introduced by Brezzi, Douglas, Marini (1985) + "Two families of mixed finite elements for Second Order Elliptic Problems" + + Following, e.g. Brezzi, Douglas, Fortin, Marini (1987) + "Efficient rectangular mixed finite elements in two and three space variables" + For rectangle K and degree j: + BDM_j(K) = [P_j(K)^2 + Span(curl(xy^{j+1}, x^{j+1}y))] x P_{j-1}(K) + + The resulting basis functions all have a curl whose polynomials are + of degree (j - 1). + + :arg deg: The element degree. + :arg dx: A tuple of sympy expressions, expanding the interval in the + x direction. Probably (1-x, x). + :arg dy: A tuple of sympy expressions, expanding the interval in the + y direction. Probably (1-y, y). + :arg x_mid: A sympy expression, probably 2*x-1. + :arg y_mid: A sympy expression, probably 2*y-1. + """ + + # For some functions, we need to multiply by a coefficient to ensure that + # the result curl is of the correct degree. The coefficient of the + # highest-order term of leg(deg, 2x-1), is binomial(2*deg, deg) coeff = binomial(2*deg, deg) / ((deg+1)*binomial(2*deg-2, deg-1)) - EL = tuple([(0, -leg(j, y_mid)*dx[0]) for j in range(deg)] + - [(coeff*-leg(deg-1, y_mid)*dy[0]*dy[1], -leg(deg, y_mid)*dx[0])] + - [(0, -leg(j, y_mid)*dx[1]) for j in range(deg)] + - [(coeff*leg(deg-1, y_mid)*dy[0]*dy[1], -leg(deg, y_mid)*dx[1])] + - [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + - [(-leg(deg, x_mid)*dy[0], coeff*-leg(deg-1, x_mid)*dx[0]*dx[1])] + - [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)] + - [(-leg(deg, x_mid)*dy[1], coeff*leg(deg-1, x_mid)*dx[0]*dx[1])]) - return EL + basis = tuple([(0, -leg(j, y_mid)*dx[0]) for j in range(deg)] + + [(coeff*-leg(deg-1, y_mid)*dy[0]*dy[1], -leg(deg, y_mid)*dx[0])] + + [(0, -leg(j, y_mid)*dx[1]) for j in range(deg)] + + [(coeff*leg(deg-1, y_mid)*dy[0]*dy[1], -leg(deg, y_mid)*dx[1])] + + [(-leg(j, x_mid)*dy[0], 0) for j in range(deg)] + + [(-leg(deg, x_mid)*dy[0], coeff*-leg(deg-1, x_mid)*dx[0]*dx[1])] + + [(-leg(j, x_mid)*dy[1], 0) for j in range(deg)] + + [(-leg(deg, x_mid)*dy[1], coeff*leg(deg-1, x_mid)*dx[0]*dx[1])]) + + return basis + +def bdmce_face_basis(deg, dx, dy, x_mid, y_mid): + """ Returns the basis functions associated with DoFs on the + faces of elements for the HCurl Brezz-Douglas-Marini element on + quadrilateral cells. -def f_lambda_1_2d(deg, dx, dy, x_mid, y_mid): - FL = [] + These were introduced by Brezzi, Douglas, Marini (1985) + "Two families of mixed finite elements for Second Order Elliptic Problems" + + Following, e.g. Brezzi, Douglas, Fortin, Marini (1987) + "Efficient rectangular mixed finite elements in two and three space variables" + For rectangle K and degree j: + BDM_j(K) = [P_j(K)^2 + Span(curl(xy^{j+1}, x^{j+1}y))] x P_{j-1}(K) + + The resulting basis functions all have a curl whose polynomials are + of degree (j - 1). + + :arg deg: The element degree. + :arg dx: A tuple of sympy expressions, expanding the interval in the + x direction. Probably (1-x, x). + :arg dy: A tuple of sympy expressions, expanding the interval in the + y direction. Probably (1-y, y). + :arg x_mid: A sympy expression, probably 2*x-1. + :arg y_mid: A sympy expression, probably 2*y-1. + """ + + basis = [] for k in range(2, deg+1): for j in range(k-1): - FL += [(0, leg(j, x_mid)*leg(k-2-j, y_mid)*dx[0]*dx[1])] - FL += [(leg(k-2-j, x_mid)*leg(j, y_mid)*dy[0]*dy[1], 0)] + basis += [(0, leg(j, x_mid)*leg(k-2-j, y_mid)*dx[0]*dx[1])] + basis += [(leg(k-2-j, x_mid)*leg(j, y_mid)*dy[0]*dy[1], 0)] - return tuple(FL) + return tuple(basis) class BrezziDouglasMariniCubeEdge(BrezziDouglasMariniCube): - def __init__(self, ref_el, degree): - if degree < 1: - raise Exception("BDMcf_k elements only valid for k >= 1") - - flat_el = flatten_reference_cube(ref_el) - dim = flat_el.get_spatial_dimension() - if dim != 2: - raise Exception("BDMcf_k elements only valid for dimension 2") + """ + The Brezzi-Douglas-Marini HCurl element on quadrilateral cells. - verts = flat_el.get_vertices() + :arg ref_el: The reference element. + :arg k: The degree. + """ - dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) - dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) - x_mid = 2*x-(verts[-1][0] + verts[0][0]) - y_mid = 2*y-(verts[-1][1] + verts[0][1]) + def __init__(self, ref_el, degree): - EL = e_lambda_1_2d(degree, dx, dy, x_mid, y_mid) - if degree >= 2: - FL = f_lambda_1_2d(degree, dx, dy, x_mid, y_mid) - else: - FL = () - bdmce_list = EL + FL + bdmce_list = construct_bdmce_basis(ref_el, degree) self.basis = {(0, 0): Array(bdmce_list)} - super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, degree=degree, mapping="covariant piola") + + super(BrezziDouglasMariniCubeEdge, self).__init__(ref_el=ref_el, degree=degree, + mapping="covariant piola") class BrezziDouglasMariniCubeFace(BrezziDouglasMariniCube): + """ + The Brezzi-Douglas-Marini HDiv element on quadrilateral cells. + + :arg ref_el: The reference element. + :arg k: The degree. + """ + def __init__(self, ref_el, degree): - if degree < 1: - raise Exception("BDMcf_k elements only valid for k >= 1") - flat_el = flatten_reference_cube(ref_el) - dim = flat_el.get_spatial_dimension() - if dim != 2: - raise Exception("BDMcf_k elements only valid for dimension 2") - - verts = flat_el.get_vertices() - - dx = ((verts[-1][0] - x)/(verts[-1][0] - verts[0][0]), (x - verts[0][0])/(verts[-1][0] - verts[0][0])) - dy = ((verts[-1][1] - y)/(verts[-1][1] - verts[0][1]), (y - verts[0][1])/(verts[-1][1] - verts[0][1])) - x_mid = 2*x-(verts[-1][0] + verts[0][0]) - y_mid = 2*y-(verts[-1][1] + verts[0][1]) - - EL = e_lambda_1_2d(degree, dx, dy, x_mid, y_mid) - if degree >= 2: - FL = f_lambda_1_2d(degree, dx, dy, x_mid, y_mid) - else: - FL = () - bdmcf_list = EL + FL - bdmcf_list = [[-a[1], a[0]] for a in bdmcf_list] + bdmce_list = construct_bdmce_basis(ref_el, degree) + + # BDMCF functions are rotations of BDMCE functions + bdmcf_list = [[-a[1], a[0]] for a in bdmce_list] self.basis = {(0, 0): Array(bdmcf_list)} - super(BrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree, mapping="contravariant piola") + super(BrezziDouglasMariniCubeFace, self).__init__(ref_el=ref_el, degree=degree, + mapping="contravariant piola") + + +def construct_bdmce_basis(ref_el, degree): + """ + Return the basis functions for a particular BDMCE space as a list. + + :arg ref_el: The reference element. + :arg k: The degree. + """ + + # Extract the vertices from the reference element + flat_el = flatten_reference_cube(ref_el) + verts = flat_el.get_vertices() + + # dx on reference quad is (1-x, x) + dx = ((verts[-1][0] - x) / (verts[-1][0] - verts[0][0]), + (x - verts[0][0]) / (verts[-1][0] - verts[0][0])) + # dy on reference quad is (1-y, y) + dy = ((verts[-1][1] - y) / (verts[-1][1] - verts[0][1]), + (y - verts[0][1]) / (verts[-1][1] - verts[0][1])) + + # x_mid and y_mid are (2x-1) and (2y-1) + x_mid = 2*x-(verts[-1][0] + verts[0][0]) + y_mid = 2*y-(verts[-1][1] + verts[0][1]) + + # Compute basis functions for BDMcE + edge_basis = bdmce_edge_basis(degree, dx, dy, x_mid, y_mid) + face_basis = bdmce_face_basis(degree, dx, dy, x_mid, y_mid) if degree > 1 else () + + return edge_basis + face_basis + + +def numpy_lambdify(X, F, modules="numpy", dummify=False): + '''Unfortunately, SymPy's own lambdify() doesn't work well with + NumPy in that simple functions like + lambda x: 1.0, + when evaluated with NumPy arrays, return just "1.0" instead of + an array of 1s with the same shape as x. This function does that. + ''' + try: + lambda_x = [numpy_lambdify(X, f, modules=modules, dummify=dummify) for f in F] + except TypeError: # 'function' object is not iterable + # SymPy's lambdify also works on functions that return arrays. + # However, use it componentwise here so we can add 0*x to each + # component individually. This is necessary to maintain shapes + # if evaluated with NumPy arrays. + lmbd_tmp = lambdify(X, F, modules=modules, dummify=dummify) + lambda_x = lambda u: lmbd_tmp(*[v for v in u]) + 0 * u[0] + + return lambda_x From c8bc3831374e688ad7afb59a3325e0517f02370f Mon Sep 17 00:00:00 2001 From: Tom Bendall Date: Thu, 28 Jan 2021 11:24:27 +0000 Subject: [PATCH 21/23] Lint fixes --- FIAT/brezzi_douglas_marini_cube.py | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index efc6878c5..5b1db8062 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -99,7 +99,6 @@ def __init__(self, ref_el, degree, mapping): self._degree = degree self.flat_el = flat_el - def degree(self): """Return the degree of the polynomial space.""" return self._degree @@ -129,12 +128,11 @@ def tabulate(self, order, points, entity=None): entity = (self.ref_el.get_dimension(), 0) entity_dim, entity_id = entity - transform = self.ref_el.get_entity_transform(entity_dim, entity_id) points = np.asarray(points) npoints, pointdim = points.shape # Turn analytic basis functions into python functions via lambdify - basis_callable = {(0, 0): numpy_lambdify(variables, self.basis[(0,0)], + basis_callable = {(0, 0): numpy_lambdify(variables, self.basis[(0, 0)], modules="numpy")} # Dictionary of values of functions @@ -146,15 +144,14 @@ def tabulate(self, order, points, entity=None): try: callable = basis_callable[alpha] except KeyError: - diff_basis = diff(self.basis[(0,0)], *zip(variables, alpha)) + diff_basis = diff(self.basis[(0, 0)], *zip(variables, alpha)) callable = numpy_lambdify(variables, diff_basis, modules="numpy") self.basis[alpha] = diff_basis # tabulate by passing points through all lambdified functions # resulting array has shape (len(self.basis), spatial_dim, npoints) T = np.array([[[func_component(point) for point in points] - for func_component in func] - for func in callable]) + for func_component in func] for func in callable]) phivals[alpha] = T @@ -172,7 +169,7 @@ def entity_closure_dofs(self): def value_shape(self): """Return the value shape of the finite element functions.""" - return np.shape(self.basis[(0,0)][0]) + return np.shape(self.basis[(0, 0)][0]) def dmats(self): raise NotImplementedError From 03edbe922ad77b5f549b5b4e07652387d88c6516 Mon Sep 17 00:00:00 2001 From: Tom Bendall Date: Thu, 28 Jan 2021 11:25:56 +0000 Subject: [PATCH 22/23] fix doc string text --- FIAT/brezzi_douglas_marini_cube.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index 5b1db8062..d4c9d4175 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -183,7 +183,7 @@ def space_dimension(self): def bdmce_edge_basis(deg, dx, dy, x_mid, y_mid): - """ Returns the basis functions associated with DoFs on the + """Returns the basis functions associated with DoFs on the edges of elements for the HCurl Brezz-Douglas-Marini element on quadrilateral cells. @@ -225,7 +225,7 @@ def bdmce_edge_basis(deg, dx, dy, x_mid, y_mid): def bdmce_face_basis(deg, dx, dy, x_mid, y_mid): - """ Returns the basis functions associated with DoFs on the + """Returns the basis functions associated with DoFs on the faces of elements for the HCurl Brezz-Douglas-Marini element on quadrilateral cells. From 5040d42bf9856b547da05d379a969f754a870f5f Mon Sep 17 00:00:00 2001 From: Tom Bendall Date: Sat, 13 Mar 2021 09:08:15 +0000 Subject: [PATCH 23/23] allow tabulation on edges --- FIAT/brezzi_douglas_marini_cube.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/FIAT/brezzi_douglas_marini_cube.py b/FIAT/brezzi_douglas_marini_cube.py index d4c9d4175..7a5ab500c 100644 --- a/FIAT/brezzi_douglas_marini_cube.py +++ b/FIAT/brezzi_douglas_marini_cube.py @@ -128,7 +128,8 @@ def tabulate(self, order, points, entity=None): entity = (self.ref_el.get_dimension(), 0) entity_dim, entity_id = entity - points = np.asarray(points) + transform = self.ref_el.get_entity_transform(entity_dim, entity_id) + points = np.asarray(list(map(transform, points))) npoints, pointdim = points.shape # Turn analytic basis functions into python functions via lambdify