From 38751784d5f1d7f3fa604a50c8cc7d181725d4ed Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Tue, 7 Feb 2017 14:56:54 +0000 Subject: [PATCH 01/12] Start adding TensorProductCell support for HDivTrace --- FIAT/hdiv_trace.py | 178 +++++++++++++++++++++++------------- test/unit/test_hdivtrace.py | 5 +- 2 files changed, 117 insertions(+), 66 deletions(-) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index ff2440538..fcded605b 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -16,14 +16,17 @@ # along with FIAT. If not, see . from __future__ import absolute_import, print_function, division +from six import iteritems import numpy as np -from FIAT.discontinuous_lagrange import DiscontinuousLagrange -from FIAT.reference_element import ufc_simplex +from FIAT.reference_element import (ufc_simplex, Point, Simplex, + TensorProductCell, + FiredrakeQuadrilateral) from FIAT.functional import PointEvaluation from FIAT.polynomial_set import mis -from FIAT import FiniteElement -from FIAT import dual_set +from FIAT import (FiniteElement, dual_set, + DiscontinuousLagrange, + TensorProductElement) # Numerical tolerance for facet-entity identifications epsilon = 1e-10 @@ -49,46 +52,54 @@ class HDivTrace(FiniteElement): """ def __init__(self, ref_el, degree): + """Constructor for the HDivTrace element. + + :arg ref_el: a reference element, which may be a tensor product + cell. + :arg degree: the degree of approximation. + """ sd = ref_el.get_spatial_dimension() if sd in (0, 1): - raise ValueError("Cannot use this trace class on a %d-dimensional cell." % sd) + raise ValueError("Cannot take the trace of a %d-dim cell." % sd) - # Constructing facet element as a discontinuous Lagrange element - dglagrange = DiscontinuousLagrange(ufc_simplex(sd - 1), degree) - - # Construct entity ids (assigning top. dim. and initializing as empty) + facet_sd = sd - 1 + dg_elements = {} entity_dofs = {} - - # Looping over dictionary of cell topology to construct the empty - # dictionary for entity ids of the trace element topology = ref_el.get_topology() - for top_dim, entities in topology.items(): + for top_dim, entities in iteritems(topology): + cell = ref_el.construct_subelement(top_dim) entity_dofs[top_dim] = {} + + if cell.get_spatial_dimension() == facet_sd: + dg_elements[top_dim] = construct_dg_element(cell, degree) + for entity in entities: entity_dofs[top_dim][entity] = [] - # Filling in entity ids and generating points for dual basis - nf = dglagrange.space_dimension() - points = [] - num_facets = sd + 1 - for f in range(num_facets): - entity_dofs[sd - 1][f] = range(f * nf, (f + 1) * nf) + offset = 0 + pts = [] + for facet_dim in dg_elements: + element = dg_elements[facet_dim] + nf = element.space_dimension() + num_facets = len(topology[facet_dim]) - for dof in dglagrange.dual_basis(): - facet_point = list(dof.get_point_dict().keys())[0] - transform = ref_el.get_entity_transform(sd - 1, f) - points.append(tuple(transform(facet_point))) + for i in range(num_facets): + entity_dofs[facet_dim][i] = range(offset + i * nf, + offset + (i + 1) * nf) + for dof in element.dual_basis(): + facet_pt = list(dof.get_point_dict().keys())[0] + transform = ref_el.get_entity_transform(facet_dim, i) + pts.append(tuple(transform(facet_pt))) - # Setting up dual basis - only point evaluations - nodes = [PointEvaluation(ref_el, pt) for pt in points] - dual = dual_set.DualSet(nodes, ref_el, entity_dofs) + offset += entity_dofs[facet_dim][num_facets - 1][-1] + 1 - super(HDivTrace, self).__init__(ref_el, dual, dglagrange.get_order(), - dglagrange.get_formdegree(), dglagrange.mapping()[0]) - # Set up facet element - self.facet_element = dglagrange + nodes = [PointEvaluation(ref_el, pt) for pt in pts] + dual = dual_set.DualSet(nodes, ref_el, entity_dofs) - # degree for quadrature rule + super(HDivTrace, self).__init__(ref_el, dual, order=degree, + formdegree=facet_sd, + mapping="affine") + self.dg_elements = dg_elements self.polydegree = degree def degree(self): @@ -125,61 +136,64 @@ def tabulate(self, order, points, entity=None): fact that performing cell-wise tabulations, or asking for any order of derivative evaluations, are not mathematically well-defined. """ - facet_dim = self.ref_el.get_spatial_dimension() - 1 - sdim = self.space_dimension() - nf = self.facet_element.space_dimension() - - # Initializing dictionary with zeros + sd = self.ref_el.get_spatial_dimension() + facet_sd = sd - 1 phivals = {} for i in range(order + 1): - alphas = mis(self.ref_el.get_spatial_dimension(), i) + alphas = mis(sd, i) + for alpha in alphas: - phivals[alpha] = np.zeros(shape=(sdim, len(points))) - evalkey = (0,) * (facet_dim + 1) + phivals[alpha] = np.zeros(shape=(self.space_dimension(), + len(points))) + + evalkey = (0,) * sd - # If entity is None, identify facet using numerical tolerance and - # return the tabulated values if entity is None: - # Attempt to identify which facet (if any) the given points are on + if not isinstance(self.ref_el, Simplex): + raise NotImplementedError( + "Tabulating this element on a %s cell without providing " + "an entity is not currently supported." % type(self.ref_el) + ) + vertices = self.ref_el.vertices coordinates = barycentric_coordinates(points, vertices) - (unique_facet, success) = extract_unique_facet(coordinates) + unique_facet, success = extract_unique_facet(coordinates) - # If successful, insert evaluations if success: - # Map points to the reference facet new_points = map_to_reference_facet(points, vertices, unique_facet) + element = self.dg_elements[facet_sd] + nf = element.space_dimension() + nonzerovals = list(element.tabulate(order, new_points).values())[0] + + phivals[evalkey][nf*unique_facet:nf * (unique_facet + 1), :] = nonzerovals - # Retrieve values by tabulating the DiscontinuousLagrange element - nonzerovals = list(self.facet_element.tabulate(order, new_points).values())[0] - phivals[evalkey][nf*unique_facet:nf*(unique_facet + 1), :] = nonzerovals - # Otherwise, return NaNs else: - for key in phivals.keys(): - phivals[key] = np.full(shape=(sdim, len(points)), fill_value=np.nan) + for key in phivals: + phivals[key] = np.full(shape=(sd, len(points)), fill_value=np.nan) return phivals entity_dim, entity_id = entity - # If the user is directly specifying cell-wise tabulation, return TraceErrors in dict for - # appropriate handling in the form compiler - if entity_dim != facet_dim: - for key in phivals.keys(): - phivals[key] = TraceError("Attempting to tabulate a %d-entity. Expecting a %d-entitiy" % (entity_dim, facet_dim)) - return phivals + if entity_dim not in self.dg_elements: + for key in phivals: + msg = "Tabulating the HDivTrace element is only allowed on facet entities" + phivals[key] = TraceError(msg) else: - # Retrieve function evaluations (order = 0 case) - nonzerovals = list(self.facet_element.tabulate(0, points).values())[0] - phivals[evalkey][nf*entity_id:nf*(entity_id + 1), :] = nonzerovals + element = self.dg_elements[entity_dim] + nf = element.space_dimension() + nonzerovals = list(element.tabulate(0, points).values())[0] + + phivals[evalkey][nf*entity_id:nf * (entity_id + 1), :] = nonzerovals - # If asking for gradient evaluations, insert TraceError in gradient evaluations if order > 0: - for key in phivals.keys(): + msg = "Gradients on trace elements are not well-defined" + for key in phivals: if key != evalkey: - phivals[key] = TraceError("Gradient evaluations are illegal on trace elements.") - return phivals + phivals[key] = TraceError(msg) + + return phivals def value_shape(self): """Return the value shape of the finite element functions.""" @@ -195,6 +209,42 @@ def get_num_members(self, arg): raise NotImplementedError("get_num_members not implemented for the trace element.") +def construct_dg_element(ref_el, degree): + """Constructs a discontinuous galerkin element of a given degree + on a particular reference cell. + """ + if isinstance(ref_el, Simplex): + dg_element = DiscontinuousLagrange(ref_el, degree) + + # Quadrilateral facets could be on a FiredrakeQuadrilateral. + # In this case, we treat this as an interval x interval cell: + elif isinstance(ref_el, FiredrakeQuadrilateral): + dg_a = DiscontinuousLagrange(ufc_simplex(1), degree) + dg_b = DiscontinuousLagrange(ufc_simplex(1), degree) + dg_element = TensorProductElement(dg_a, dg_b) + + # This handles the more general case for facets: + elif isinstance(ref_el, TensorProductCell): + A, B = ref_el.cells + + if isinstance(A, Point) and not isinstance(B, Point): + dg_element = construct_dg_element(B, degree) + + elif isinstance(B, Point) and not isinstance(A, Point): + dg_element = construct_dg_element(A, degree) + + else: + dg_a = construct_dg_element(A, degree) + dg_b = construct_dg_element(B, degree) + dg_element = TensorProductElement(dg_a, dg_b) + else: + raise NotImplementedError( + "Reference cells of type %s not currently supported" % type(ref_el) + ) + + return dg_element + + # The following functions are credited to Marie E. Rognes: def extract_unique_facet(coordinates, tolerance=epsilon): """Determines whether a set of points (described in its barycentric coordinates) diff --git a/test/unit/test_hdivtrace.py b/test/unit/test_hdivtrace.py index cfd96d976..5138fa520 100644 --- a/test/unit/test_hdivtrace.py +++ b/test/unit/test_hdivtrace.py @@ -43,7 +43,8 @@ def test_basis_values(dim, degree): ref_el = ufc_simplex(dim) quadrule = make_quadrature(ufc_simplex(dim - 1), degree + 1) fiat_element = HDivTrace(ref_el, degree) - nf = fiat_element.facet_element.space_dimension() + facet_element = fiat_element.dg_elements[dim - 1] + nf = facet_element.space_dimension() for facet_id in range(dim + 1): # Tabulate without an entity pair given --- need to map to cell coordinates @@ -58,7 +59,7 @@ def test_basis_values(dim, degree): for test_degree in range(degree + 1): coeffs = [n(lambda x: x[0]**test_degree) - for n in fiat_element.facet_element.dual.nodes] + for n in facet_element.dual.nodes] cintegral = np.dot(coeffs, np.dot(ctab, quadrule.wts)) eintegral = np.dot(coeffs, np.dot(etab, quadrule.wts)) From d32c67f68f976494126227352906f1da35ef7fcf Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Tue, 7 Feb 2017 17:01:57 +0000 Subject: [PATCH 02/12] Test the trace on a quad --- test/unit/test_hdivtrace.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/test/unit/test_hdivtrace.py b/test/unit/test_hdivtrace.py index 5138fa520..4b0dcf476 100644 --- a/test/unit/test_hdivtrace.py +++ b/test/unit/test_hdivtrace.py @@ -71,6 +71,36 @@ def test_basis_values(dim, degree): assert np.allclose(eintegral, reference, rtol=1e-14) +@pytest.mark.parametrize("degree", range(4)) +def test_quad_trace(degree): + """Test the trace element defined on a quadrilateral cell""" + from FIAT import ufc_simplex, HDivTrace, make_quadrature + from FIAT.reference_element import TensorProductCell + + tpc = TensorProductCell(ufc_simplex(1), ufc_simplex(1)) + fiat_element = HDivTrace(tpc, degree) + facet_elements = fiat_element.dg_elements + quadrule = make_quadrature(ufc_simplex(1), degree + 1) + + for entity in [((0, 1), 0), ((0, 1), 1), ((1, 0), 0), ((1, 0), 1)]: + entity_dim, entity_id = entity + element = facet_elements[entity_dim] + nf = element.space_dimension() + + tab = fiat_element.tabulate(0, quadrule.pts, + entity)[(0, 0)][nf*entity_id:nf*(entity_id + 1)] + + for test_degree in range(degree + 1): + coeffs = [n(lambda x: x[0]**test_degree) + for n in element.dual.nodes] + + integral = np.dot(coeffs, np.dot(tab, quadrule.wts)) + + reference = np.dot([x[0]**test_degree + for x in quadrule.pts], quadrule.wts) + assert np.allclose(integral, reference, rtol=1e-14) + + @pytest.mark.parametrize("dim", (2, 3)) @pytest.mark.parametrize("order", range(1, 4)) @pytest.mark.parametrize("degree", range(4)) From 68bc9ee283ea07387a79587f091f80cda6f7e663 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 8 Feb 2017 16:12:04 +0000 Subject: [PATCH 03/12] Add multi-degree support for tensor product cells --- FIAT/hdiv_trace.py | 68 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 50 insertions(+), 18 deletions(-) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index fcded605b..4d5023fbd 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -54,14 +54,28 @@ class HDivTrace(FiniteElement): def __init__(self, ref_el, degree): """Constructor for the HDivTrace element. - :arg ref_el: a reference element, which may be a tensor product + :arg ref_el: A reference element, which may be a tensor product cell. - :arg degree: the degree of approximation. + :arg degree: The degree of approximation. If on a tensor product + cell, then provide a tuple of degrees. """ sd = ref_el.get_spatial_dimension() if sd in (0, 1): raise ValueError("Cannot take the trace of a %d-dim cell." % sd) + if isinstance(ref_el, TensorProductCell): + if isinstance(degree, tuple): + spanning_degrees = degree + + else: + spanning_degrees = (degree,) * len(ref_el.cells) + else: + assert isinstance(ref_el, (Simplex, FiredrakeQuadrilateral)) + assert not isinstance(degree, tuple), ( + "Must have a tensor product cell if providing multiple degrees" + ) + spanning_degrees = (degree,) + facet_sd = sd - 1 dg_elements = {} entity_dofs = {} @@ -71,14 +85,15 @@ def __init__(self, ref_el, degree): entity_dofs[top_dim] = {} if cell.get_spatial_dimension() == facet_sd: - dg_elements[top_dim] = construct_dg_element(cell, degree) + dg_elements[top_dim] = construct_dg_element(cell, + spanning_degrees) for entity in entities: entity_dofs[top_dim][entity] = [] offset = 0 pts = [] - for facet_dim in dg_elements: + for facet_dim in sorted(dg_elements): element = dg_elements[facet_dim] nf = element.space_dimension() num_facets = len(topology[facet_dim]) @@ -96,11 +111,14 @@ def __init__(self, ref_el, degree): nodes = [PointEvaluation(ref_el, pt) for pt in pts] dual = dual_set.DualSet(nodes, ref_el, entity_dofs) - super(HDivTrace, self).__init__(ref_el, dual, order=degree, + deg = max([e.degree() for e in dg_elements.values()]) + + super(HDivTrace, self).__init__(ref_el, dual, order=deg, formdegree=facet_sd, mapping="affine") self.dg_elements = dg_elements - self.polydegree = degree + self.polydegree = deg + self._spanning_degrees = spanning_degrees def degree(self): """Return the degree of the (embedding) polynomial space.""" @@ -160,16 +178,22 @@ def tabulate(self, order, points, entity=None): unique_facet, success = extract_unique_facet(coordinates) if success: - new_points = map_to_reference_facet(points, vertices, unique_facet) + new_points = map_to_reference_facet(points, + vertices, + unique_facet) element = self.dg_elements[facet_sd] nf = element.space_dimension() - nonzerovals = list(element.tabulate(order, new_points).values())[0] + nonzerovals = list(element.tabulate(order, + new_points).values())[0] - phivals[evalkey][nf*unique_facet:nf * (unique_facet + 1), :] = nonzerovals + sindex = nf * unique_facet + eindex = nf * (unique_facet + 1) + phivals[evalkey][sindex:eindex, :] = nonzerovals else: for key in phivals: - phivals[key] = np.full(shape=(sd, len(points)), fill_value=np.nan) + phivals[key] = np.full(shape=(sd, len(points)), + fill_value=np.nan) return phivals @@ -177,7 +201,7 @@ def tabulate(self, order, points, entity=None): if entity_dim not in self.dg_elements: for key in phivals: - msg = "Tabulating the HDivTrace element is only allowed on facet entities" + msg = "The HDivTrace element can only be tabulated on facets." phivals[key] = TraceError(msg) else: @@ -185,10 +209,12 @@ def tabulate(self, order, points, entity=None): nf = element.space_dimension() nonzerovals = list(element.tabulate(0, points).values())[0] - phivals[evalkey][nf*entity_id:nf * (entity_id + 1), :] = nonzerovals + sindex = nf * entity_id + eindex = nf * (entity_id + 1) + phivals[evalkey][sindex:eindex, :] = nonzerovals if order > 0: - msg = "Gradients on trace elements are not well-defined" + msg = "Gradients on trace elements are not well-defined." for key in phivals: if key != evalkey: phivals[key] = TraceError(msg) @@ -209,33 +235,39 @@ def get_num_members(self, arg): raise NotImplementedError("get_num_members not implemented for the trace element.") -def construct_dg_element(ref_el, degree): +def construct_dg_element(ref_el, spanning_degrees): """Constructs a discontinuous galerkin element of a given degree on a particular reference cell. """ if isinstance(ref_el, Simplex): + degree, = spanning_degrees dg_element = DiscontinuousLagrange(ref_el, degree) # Quadrilateral facets could be on a FiredrakeQuadrilateral. # In this case, we treat this as an interval x interval cell: elif isinstance(ref_el, FiredrakeQuadrilateral): + degree, = spanning_degrees dg_a = DiscontinuousLagrange(ufc_simplex(1), degree) dg_b = DiscontinuousLagrange(ufc_simplex(1), degree) dg_element = TensorProductElement(dg_a, dg_b) # This handles the more general case for facets: elif isinstance(ref_el, TensorProductCell): + assert len(spanning_degrees) == 2, ( + "Must provide two degrees for tensor product cells" + ) + d1, d2 = spanning_degrees A, B = ref_el.cells if isinstance(A, Point) and not isinstance(B, Point): - dg_element = construct_dg_element(B, degree) + dg_element = construct_dg_element(B, (d2,)) elif isinstance(B, Point) and not isinstance(A, Point): - dg_element = construct_dg_element(A, degree) + dg_element = construct_dg_element(A, (d1,)) else: - dg_a = construct_dg_element(A, degree) - dg_b = construct_dg_element(B, degree) + dg_a = construct_dg_element(A, (d1,)) + dg_b = construct_dg_element(B, (d2,)) dg_element = TensorProductElement(dg_a, dg_b) else: raise NotImplementedError( From 2e0fb6dc59bc0c2612b83e5aa53d9d06fee88d90 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 9 Feb 2017 12:31:32 +0000 Subject: [PATCH 04/12] Fix offset in Trace tabulation --- FIAT/hdiv_trace.py | 73 +++++++++++++++++++++++++++++++++++++++------- 1 file changed, 62 insertions(+), 11 deletions(-) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index 4d5023fbd..f28644042 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -16,7 +16,7 @@ # along with FIAT. If not, see . from __future__ import absolute_import, print_function, division -from six import iteritems +from six import iteritems, itervalues import numpy as np from FIAT.reference_element import (ufc_simplex, Point, Simplex, @@ -57,25 +57,36 @@ def __init__(self, ref_el, degree): :arg ref_el: A reference element, which may be a tensor product cell. :arg degree: The degree of approximation. If on a tensor product - cell, then provide a tuple of degrees. + cell, then provide a tuple of degrees if you want + varying degrees. """ sd = ref_el.get_spatial_dimension() if sd in (0, 1): raise ValueError("Cannot take the trace of a %d-dim cell." % sd) + # Store the spanning degrees if on a tensor product cell if isinstance(ref_el, TensorProductCell): if isinstance(degree, tuple): + # It is possible to have varying degrees in the component + # elements of a tensor product element. This is the same + # for the trace element defined on a tensor product cell. spanning_degrees = degree else: + # If one degree is provided, we splat the degree + # over all component cells spanning_degrees = (degree,) * len(ref_el.cells) else: assert isinstance(ref_el, (Simplex, FiredrakeQuadrilateral)) + + # Cannot have varying degrees for these reference cells assert not isinstance(degree, tuple), ( "Must have a tensor product cell if providing multiple degrees" ) spanning_degrees = (degree,) + # Initialize entity dofs and construct the DG elements + # for the facets facet_sd = sd - 1 dg_elements = {} entity_dofs = {} @@ -84,13 +95,16 @@ def __init__(self, ref_el, degree): cell = ref_el.construct_subelement(top_dim) entity_dofs[top_dim] = {} + # We have a facet entity! if cell.get_spatial_dimension() == facet_sd: dg_elements[top_dim] = construct_dg_element(cell, spanning_degrees) - + # Initialize for entity in entities: entity_dofs[top_dim][entity] = [] + # Compute the dof numbering for all facet entities + # and extract nodes offset = 0 pts = [] for facet_dim in sorted(dg_elements): @@ -101,6 +115,8 @@ def __init__(self, ref_el, degree): for i in range(num_facets): entity_dofs[facet_dim][i] = range(offset + i * nf, offset + (i + 1) * nf) + + # Run over nodes and collect the points for point evaluations for dof in element.dual_basis(): facet_pt = list(dof.get_point_dict().keys())[0] transform = ref_el.get_entity_transform(facet_dim, i) @@ -108,16 +124,24 @@ def __init__(self, ref_el, degree): offset += entity_dofs[facet_dim][num_facets - 1][-1] + 1 + # Setting up dual basis - only point evaluations nodes = [PointEvaluation(ref_el, pt) for pt in pts] dual = dual_set.DualSet(nodes, ref_el, entity_dofs) + # Degree of the element deg = max([e.degree() for e in dg_elements.values()]) super(HDivTrace, self).__init__(ref_el, dual, order=deg, formdegree=facet_sd, mapping="affine") + + # Set up facet elements self.dg_elements = dg_elements + + # Degree for quadrature rule self.polydegree = deg + + # Store the spanning degrees self._spanning_degrees = spanning_degrees def degree(self): @@ -156,6 +180,8 @@ def tabulate(self, order, points, entity=None): """ sd = self.ref_el.get_spatial_dimension() facet_sd = sd - 1 + + # Initializing dictionary with zeros phivals = {} for i in range(order + 1): alphas = mis(sd, i) @@ -166,30 +192,39 @@ def tabulate(self, order, points, entity=None): evalkey = (0,) * sd + # If entity is None, identify facet using numerical tolerance and + # return the tabulated values if entity is None: + # NOTE: Numerical approximation of the facet id is currently only + # implemented for simplex reference cells. if not isinstance(self.ref_el, Simplex): raise NotImplementedError( "Tabulating this element on a %s cell without providing " "an entity is not currently supported." % type(self.ref_el) ) + # Attempt to identify which facet (if any) the given points are on vertices = self.ref_el.vertices coordinates = barycentric_coordinates(points, vertices) unique_facet, success = extract_unique_facet(coordinates) + # If successful, insert evaluations if success: + # Map points to the reference facet new_points = map_to_reference_facet(points, vertices, unique_facet) + + # Retrieve values by tabulating the DG element element = self.dg_elements[facet_sd] nf = element.space_dimension() - nonzerovals = list(element.tabulate(order, - new_points).values())[0] + nonzerovals, = itervalues(element.tabulate(order, new_points)) sindex = nf * unique_facet eindex = nf * (unique_facet + 1) phivals[evalkey][sindex:eindex, :] = nonzerovals + # Otherwise, return NaNs else: for key in phivals: phivals[key] = np.full(shape=(sd, len(points)), @@ -199,20 +234,36 @@ def tabulate(self, order, points, entity=None): entity_dim, entity_id = entity + # If the user is directly specifying cell-wise tabulation, return + # TraceErrors in dict for appropriate handling in the form compiler if entity_dim not in self.dg_elements: for key in phivals: msg = "The HDivTrace element can only be tabulated on facets." phivals[key] = TraceError(msg) else: - element = self.dg_elements[entity_dim] - nf = element.space_dimension() - nonzerovals = list(element.tabulate(0, points).values())[0] + # Retrieve function evaluations (order = 0 case) + # NOTE: Depending on the cell, we need to off set + # the facet index when tabulating particular facet + # entities + offset = 0 + for facet_dim in sorted(self.dg_elements): + element = self.dg_elements[facet_dim] + nf = element.space_dimension() + num_facets = len(self.ref_el.get_topology()[facet_dim]) + + # Loop over the number of facets until we find a facet + # with matching dimension and id + for i in range(num_facets): + # Found it! + if facet_dim == entity_dim and i == entity_id: + nonzerovals, = itervalues(element.tabulate(0, points)) + phivals[evalkey][offset:offset+nf, :] = nonzerovals - sindex = nf * entity_id - eindex = nf * (entity_id + 1) - phivals[evalkey][sindex:eindex, :] = nonzerovals + offset += nf + # If asking for gradient evaluations, insert TraceError in + # gradient slots if order > 0: msg = "Gradients on trace elements are not well-defined." for key in phivals: From c6fb3c8c1de92fd7be0f5c1a6bc32f55b1e2647d Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 9 Feb 2017 12:31:56 +0000 Subject: [PATCH 05/12] Update test for HDivTrace --- test/unit/test_hdivtrace.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/unit/test_hdivtrace.py b/test/unit/test_hdivtrace.py index 4b0dcf476..71c8e873a 100644 --- a/test/unit/test_hdivtrace.py +++ b/test/unit/test_hdivtrace.py @@ -78,17 +78,18 @@ def test_quad_trace(degree): from FIAT.reference_element import TensorProductCell tpc = TensorProductCell(ufc_simplex(1), ufc_simplex(1)) - fiat_element = HDivTrace(tpc, degree) + fiat_element = HDivTrace(tpc, (degree, degree)) facet_elements = fiat_element.dg_elements quadrule = make_quadrature(ufc_simplex(1), degree + 1) - for entity in [((0, 1), 0), ((0, 1), 1), ((1, 0), 0), ((1, 0), 1)]: - entity_dim, entity_id = entity + for i, entity in enumerate([((0, 1), 0), ((0, 1), 1), + ((1, 0), 0), ((1, 0), 1)]): + entity_dim, _ = entity element = facet_elements[entity_dim] nf = element.space_dimension() tab = fiat_element.tabulate(0, quadrule.pts, - entity)[(0, 0)][nf*entity_id:nf*(entity_id + 1)] + entity)[(0, 0)][nf*i:nf*(i+1)] for test_degree in range(degree + 1): coeffs = [n(lambda x: x[0]**test_degree) From de9bfae6eeff1379f418ccb098b4b2114a332683 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 29 Mar 2017 13:54:53 +0100 Subject: [PATCH 06/12] Add a more robust pattern for trace degrees --- FIAT/hdiv_trace.py | 52 +++++++++++++++++++++------------------------- 1 file changed, 24 insertions(+), 28 deletions(-) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index f28644042..1daf6179a 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -21,7 +21,8 @@ import numpy as np from FIAT.reference_element import (ufc_simplex, Point, Simplex, TensorProductCell, - FiredrakeQuadrilateral) + FiredrakeQuadrilateral, + TENSORPRODUCT) from FIAT.functional import PointEvaluation from FIAT.polynomial_set import mis from FIAT import (FiniteElement, dual_set, @@ -64,18 +65,17 @@ def __init__(self, ref_el, degree): if sd in (0, 1): raise ValueError("Cannot take the trace of a %d-dim cell." % sd) - # Store the spanning degrees if on a tensor product cell - if isinstance(ref_el, TensorProductCell): - if isinstance(degree, tuple): - # It is possible to have varying degrees in the component - # elements of a tensor product element. This is the same - # for the trace element defined on a tensor product cell. - spanning_degrees = degree + # Store the degrees if on a tensor product cell + if ref_el.get_shape() == TENSORPRODUCT: + try: + degree = tuple(degree) + except TypeError: + degree = (degree,) * len(ref_el.cells) - else: - # If one degree is provided, we splat the degree - # over all component cells - spanning_degrees = (degree,) * len(ref_el.cells) + assert len(ref_el.cells) == len(degree), ( + "Number of specified degrees must be equal " + "to the number of cells." + ) else: assert isinstance(ref_el, (Simplex, FiredrakeQuadrilateral)) @@ -83,7 +83,7 @@ def __init__(self, ref_el, degree): assert not isinstance(degree, tuple), ( "Must have a tensor product cell if providing multiple degrees" ) - spanning_degrees = (degree,) + degree = (degree,) # Initialize entity dofs and construct the DG elements # for the facets @@ -97,8 +97,7 @@ def __init__(self, ref_el, degree): # We have a facet entity! if cell.get_spatial_dimension() == facet_sd: - dg_elements[top_dim] = construct_dg_element(cell, - spanning_degrees) + dg_elements[top_dim] = construct_dg_element(cell, degree) # Initialize for entity in entities: entity_dofs[top_dim][entity] = [] @@ -141,9 +140,6 @@ def __init__(self, ref_el, degree): # Degree for quadrature rule self.polydegree = deg - # Store the spanning degrees - self._spanning_degrees = spanning_degrees - def degree(self): """Return the degree of the (embedding) polynomial space.""" return self.polydegree @@ -172,11 +168,11 @@ def tabulate(self, order, points, entity=None): .. note :: - Performing illegal tabulations on this element will result in either - a tabulation table of `numpy.nan` arrays (`entity=None` case), or - insertions of the `TraceError` exception class. This is due to the - fact that performing cell-wise tabulations, or asking for any order - of derivative evaluations, are not mathematically well-defined. + Performing illegal tabulations on this element will result in either + a tabulation table of `numpy.nan` arrays (`entity=None` case), or + insertions of the `TraceError` exception class. This is due to the + fact that performing cell-wise tabulations, or asking for any order + of derivative evaluations, are not mathematically well-defined. """ sd = self.ref_el.get_spatial_dimension() facet_sd = sd - 1 @@ -286,28 +282,28 @@ def get_num_members(self, arg): raise NotImplementedError("get_num_members not implemented for the trace element.") -def construct_dg_element(ref_el, spanning_degrees): +def construct_dg_element(ref_el, degrees): """Constructs a discontinuous galerkin element of a given degree on a particular reference cell. """ if isinstance(ref_el, Simplex): - degree, = spanning_degrees + degree, = degrees dg_element = DiscontinuousLagrange(ref_el, degree) # Quadrilateral facets could be on a FiredrakeQuadrilateral. # In this case, we treat this as an interval x interval cell: elif isinstance(ref_el, FiredrakeQuadrilateral): - degree, = spanning_degrees + degree, = degrees dg_a = DiscontinuousLagrange(ufc_simplex(1), degree) dg_b = DiscontinuousLagrange(ufc_simplex(1), degree) dg_element = TensorProductElement(dg_a, dg_b) # This handles the more general case for facets: elif isinstance(ref_el, TensorProductCell): - assert len(spanning_degrees) == 2, ( + assert len(degrees) == 2, ( "Must provide two degrees for tensor product cells" ) - d1, d2 = spanning_degrees + d1, d2 = degrees A, B = ref_el.cells if isinstance(A, Point) and not isinstance(B, Point): From 89919e13d8f14b8c0370d349994becc17ab262fd Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 29 Mar 2017 14:28:43 +0100 Subject: [PATCH 07/12] Simplify logic for constructing facet elements --- FIAT/hdiv_trace.py | 62 +++++++++++++++++++++++----------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index 1daf6179a..7b1b7e158 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -17,11 +17,12 @@ from __future__ import absolute_import, print_function, division from six import iteritems, itervalues +from six.moves import range import numpy as np -from FIAT.reference_element import (ufc_simplex, Point, Simplex, - TensorProductCell, - FiredrakeQuadrilateral, +from FIAT.reference_element import (ufc_simplex, POINT, + LINE, QUADRILATERAL, + TRIANGLE, TETRAHEDRON, TENSORPRODUCT) from FIAT.functional import PointEvaluation from FIAT.polynomial_set import mis @@ -77,13 +78,17 @@ def __init__(self, ref_el, degree): "to the number of cells." ) else: - assert isinstance(ref_el, (Simplex, FiredrakeQuadrilateral)) - + if not ref_el.get_shape() in [TRIANGLE, TETRAHEDRON, + QUADRILATERAL]: + raise NotImplementedError( + "Trace element on a %s not implemented" % type(ref_el) + ) # Cannot have varying degrees for these reference cells - assert not isinstance(degree, tuple), ( - "Must have a tensor product cell if providing multiple degrees" - ) - degree = (degree,) + if isinstance(degree, tuple): + raise ValueError( + "Must have a tensor product cell " + "if providing multiple degrees" + ) # Initialize entity dofs and construct the DG elements # for the facets @@ -117,7 +122,7 @@ def __init__(self, ref_el, degree): # Run over nodes and collect the points for point evaluations for dof in element.dual_basis(): - facet_pt = list(dof.get_point_dict().keys())[0] + facet_pt, = dof.get_point_dict() transform = ref_el.get_entity_transform(facet_dim, i) pts.append(tuple(transform(facet_pt))) @@ -193,7 +198,7 @@ def tabulate(self, order, points, entity=None): if entity is None: # NOTE: Numerical approximation of the facet id is currently only # implemented for simplex reference cells. - if not isinstance(self.ref_el, Simplex): + if not self.ref_el.get_shape() in [TRIANGLE, TETRAHEDRON]: raise NotImplementedError( "Tabulating this element on a %s cell without providing " "an entity is not currently supported." % type(self.ref_el) @@ -270,7 +275,7 @@ def tabulate(self, order, points, entity=None): def value_shape(self): """Return the value shape of the finite element functions.""" - return self.facet_element.value_shape() + return () def dmats(self): """Return dmats: expansion coefficients for basis function @@ -282,40 +287,35 @@ def get_num_members(self, arg): raise NotImplementedError("get_num_members not implemented for the trace element.") -def construct_dg_element(ref_el, degrees): +def construct_dg_element(ref_el, degree): """Constructs a discontinuous galerkin element of a given degree on a particular reference cell. """ - if isinstance(ref_el, Simplex): - degree, = degrees + if ref_el.get_shape() in [LINE, TRIANGLE]: dg_element = DiscontinuousLagrange(ref_el, degree) # Quadrilateral facets could be on a FiredrakeQuadrilateral. # In this case, we treat this as an interval x interval cell: - elif isinstance(ref_el, FiredrakeQuadrilateral): - degree, = degrees + elif ref_el.get_shape() == QUADRILATERAL: dg_a = DiscontinuousLagrange(ufc_simplex(1), degree) dg_b = DiscontinuousLagrange(ufc_simplex(1), degree) dg_element = TensorProductElement(dg_a, dg_b) # This handles the more general case for facets: - elif isinstance(ref_el, TensorProductCell): - assert len(degrees) == 2, ( - "Must provide two degrees for tensor product cells" + elif ref_el.get_shape() == TENSORPRODUCT: + assert len(degree) == len(ref_el.cells), ( + "Must provide the same number of degrees as the number " + "of cells that make up the tensor product cell." ) - d1, d2 = degrees - A, B = ref_el.cells - - if isinstance(A, Point) and not isinstance(B, Point): - dg_element = construct_dg_element(B, (d2,)) - - elif isinstance(B, Point) and not isinstance(A, Point): - dg_element = construct_dg_element(A, (d1,)) + sub_elements = [construct_dg_element(c, d) + for c, d in zip(ref_el.cells, degree) + if c.get_shape() != POINT] + if len(sub_elements) > 1: + dg_element = TensorProductElement(*sub_elements) else: - dg_a = construct_dg_element(A, (d1,)) - dg_b = construct_dg_element(B, (d2,)) - dg_element = TensorProductElement(dg_a, dg_b) + dg_element, = sub_elements + else: raise NotImplementedError( "Reference cells of type %s not currently supported" % type(ref_el) From 85ff2074913cfdc70f421a9975948af0f3fc94b3 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 29 Mar 2017 18:16:35 +0100 Subject: [PATCH 08/12] Fix entity dof labeling --- FIAT/hdiv_trace.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index 7b1b7e158..ec8272df1 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -117,8 +117,8 @@ def __init__(self, ref_el, degree): num_facets = len(topology[facet_dim]) for i in range(num_facets): - entity_dofs[facet_dim][i] = range(offset + i * nf, - offset + (i + 1) * nf) + entity_dofs[facet_dim][i] = list(range(offset + i * nf, + offset + (i + 1) * nf)) # Run over nodes and collect the points for point evaluations for dof in element.dual_basis(): @@ -126,7 +126,7 @@ def __init__(self, ref_el, degree): transform = ref_el.get_entity_transform(facet_dim, i) pts.append(tuple(transform(facet_pt))) - offset += entity_dofs[facet_dim][num_facets - 1][-1] + 1 + offset += nf * num_facets # Setting up dual basis - only point evaluations nodes = [PointEvaluation(ref_el, pt) for pt in pts] From ba5fa686342cc2f75ac2c2bd02d908bc1d43f480 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Tue, 4 Apr 2017 21:51:13 +0100 Subject: [PATCH 09/12] Clean up line-breaks and organize imports --- FIAT/hdiv_trace.py | 36 ++++++++++++++---------------------- 1 file changed, 14 insertions(+), 22 deletions(-) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index ec8272df1..1b55cc7b7 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -20,15 +20,16 @@ from six.moves import range import numpy as np +from FIAT.discontinuous_lagrange import DiscontinuousLagrange +from FIAT.dual_set import DualSet +from FIAT.finite_element import FiniteElement +from FIAT.functional import PointEvaluation +from FIAT.polynomial_set import mis from FIAT.reference_element import (ufc_simplex, POINT, LINE, QUADRILATERAL, TRIANGLE, TETRAHEDRON, TENSORPRODUCT) -from FIAT.functional import PointEvaluation -from FIAT.polynomial_set import mis -from FIAT import (FiniteElement, dual_set, - DiscontinuousLagrange, - TensorProductElement) +from FIAT.tensor_product import TensorProductElement # Numerical tolerance for facet-entity identifications epsilon = 1e-10 @@ -74,21 +75,16 @@ def __init__(self, ref_el, degree): degree = (degree,) * len(ref_el.cells) assert len(ref_el.cells) == len(degree), ( - "Number of specified degrees must be equal " - "to the number of cells." + "Number of specified degrees must be equal to the number of cells." ) else: - if not ref_el.get_shape() in [TRIANGLE, TETRAHEDRON, - QUADRILATERAL]: + if ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON, QUADRILATERAL]: raise NotImplementedError( "Trace element on a %s not implemented" % type(ref_el) ) # Cannot have varying degrees for these reference cells if isinstance(degree, tuple): - raise ValueError( - "Must have a tensor product cell " - "if providing multiple degrees" - ) + raise ValueError("Must have a tensor product cell if providing multiple degrees") # Initialize entity dofs and construct the DG elements # for the facets @@ -130,7 +126,7 @@ def __init__(self, ref_el, degree): # Setting up dual basis - only point evaluations nodes = [PointEvaluation(ref_el, pt) for pt in pts] - dual = dual_set.DualSet(nodes, ref_el, entity_dofs) + dual = DualSet(nodes, ref_el, entity_dofs) # Degree of the element deg = max([e.degree() for e in dg_elements.values()]) @@ -188,8 +184,7 @@ def tabulate(self, order, points, entity=None): alphas = mis(sd, i) for alpha in alphas: - phivals[alpha] = np.zeros(shape=(self.space_dimension(), - len(points))) + phivals[alpha] = np.zeros(shape=(self.space_dimension(), len(points))) evalkey = (0,) * sd @@ -198,7 +193,7 @@ def tabulate(self, order, points, entity=None): if entity is None: # NOTE: Numerical approximation of the facet id is currently only # implemented for simplex reference cells. - if not self.ref_el.get_shape() in [TRIANGLE, TETRAHEDRON]: + if self.ref_el.get_shape() not in [TRIANGLE, TETRAHEDRON]: raise NotImplementedError( "Tabulating this element on a %s cell without providing " "an entity is not currently supported." % type(self.ref_el) @@ -212,9 +207,7 @@ def tabulate(self, order, points, entity=None): # If successful, insert evaluations if success: # Map points to the reference facet - new_points = map_to_reference_facet(points, - vertices, - unique_facet) + new_points = map_to_reference_facet(points, vertices, unique_facet) # Retrieve values by tabulating the DG element element = self.dg_elements[facet_sd] @@ -228,8 +221,7 @@ def tabulate(self, order, points, entity=None): # Otherwise, return NaNs else: for key in phivals: - phivals[key] = np.full(shape=(sd, len(points)), - fill_value=np.nan) + phivals[key] = np.full(shape=(sd, len(points)), fill_value=np.nan) return phivals From 97cc5e4fe083c6e180c988f643ea3c9ef560835e Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Wed, 5 Apr 2017 10:59:14 +0100 Subject: [PATCH 10/12] Update next.rst --- doc/sphinx/source/releases/next.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/doc/sphinx/source/releases/next.rst b/doc/sphinx/source/releases/next.rst index cdcd4a352..aaa3b6749 100644 --- a/doc/sphinx/source/releases/next.rst +++ b/doc/sphinx/source/releases/next.rst @@ -11,6 +11,13 @@ Summary of changes be published (and renamed) to list the most important changes in the new release. +- Extended the discontinuous trace element ``HDivTrace`` to support tensor + product reference cells. Tabulating the trace defined on a tensor product + cell relies on the argument ``entity`` to specify a facet of the cell. The + backwards compatibility case ``entity=None`` does not support tensor product + tabulation as a result. Tabulating the trace of triangles or tetrahedron + remains unaffected and works as usual with or without an entity argument. + Detailed changes ================ From 59ccfd3d806517f5c79d37bd73df262e25a917e0 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 6 Apr 2017 12:16:35 +0100 Subject: [PATCH 11/12] Apply minor logic edits --- FIAT/hdiv_trace.py | 95 +++++++++++++++++++++++----------------------- 1 file changed, 47 insertions(+), 48 deletions(-) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index 1b55cc7b7..8e49fea75 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -113,8 +113,8 @@ def __init__(self, ref_el, degree): num_facets = len(topology[facet_dim]) for i in range(num_facets): - entity_dofs[facet_dim][i] = list(range(offset + i * nf, - offset + (i + 1) * nf)) + entity_dofs[facet_dim][i] = list(range(offset, offset + nf)) + offset += nf # Run over nodes and collect the points for point evaluations for dof in element.dual_basis(): @@ -122,8 +122,6 @@ def __init__(self, ref_el, degree): transform = ref_el.get_entity_transform(facet_dim, i) pts.append(tuple(transform(facet_pt))) - offset += nf * num_facets - # Setting up dual basis - only point evaluations nodes = [PointEvaluation(ref_el, pt) for pt in pts] dual = DualSet(nodes, ref_el, entity_dofs) @@ -182,7 +180,6 @@ def tabulate(self, order, points, entity=None): phivals = {} for i in range(order + 1): alphas = mis(sd, i) - for alpha in alphas: phivals[alpha] = np.zeros(shape=(self.space_dimension(), len(points))) @@ -204,8 +201,15 @@ def tabulate(self, order, points, entity=None): coordinates = barycentric_coordinates(points, vertices) unique_facet, success = extract_unique_facet(coordinates) - # If successful, insert evaluations - if success: + # If not successful, return NaNs + if not success: + for key in phivals: + phivals[key] = np.full(shape=(sd, len(points)), fill_value=np.nan) + + return phivals + + # Otherwise, extract non-zero values and insertion indices + else: # Map points to the reference facet new_points = map_to_reference_facet(points, vertices, unique_facet) @@ -213,55 +217,50 @@ def tabulate(self, order, points, entity=None): element = self.dg_elements[facet_sd] nf = element.space_dimension() nonzerovals, = itervalues(element.tabulate(order, new_points)) - sindex = nf * unique_facet eindex = nf * (unique_facet + 1) - phivals[evalkey][sindex:eindex, :] = nonzerovals - # Otherwise, return NaNs - else: - for key in phivals: - phivals[key] = np.full(shape=(sd, len(points)), fill_value=np.nan) + else: + entity_dim, _ = entity - return phivals + # If the user is directly specifying cell-wise tabulation, return + # TraceErrors in dict for appropriate handling in the form compiler + if entity_dim not in self.dg_elements: + for key in phivals: + msg = "The HDivTrace element can only be tabulated on facets." + phivals[key] = TraceError(msg) - entity_dim, entity_id = entity + return phivals - # If the user is directly specifying cell-wise tabulation, return - # TraceErrors in dict for appropriate handling in the form compiler - if entity_dim not in self.dg_elements: + else: + # Retrieve function evaluations (order = 0 case) + offset = 0 + for facet_dim in sorted(self.dg_elements): + element = self.dg_elements[facet_dim] + nf = element.space_dimension() + num_facets = len(self.ref_el.get_topology()[facet_dim]) + + # Loop over the number of facets until we find a facet + # with matching dimension and id + for i in range(num_facets): + # Found it! Grab insertion indices + if (facet_dim, i) == entity: + nonzerovals, = itervalues(element.tabulate(0, points)) + sindex = offset + eindex = offset + nf + + offset += nf + + # If asking for gradient evaluations, insert TraceError in + # gradient slots + if order > 0: + msg = "Gradients on trace elements are not well-defined." for key in phivals: - msg = "The HDivTrace element can only be tabulated on facets." - phivals[key] = TraceError(msg) + if key != evalkey: + phivals[key] = TraceError(msg) - else: - # Retrieve function evaluations (order = 0 case) - # NOTE: Depending on the cell, we need to off set - # the facet index when tabulating particular facet - # entities - offset = 0 - for facet_dim in sorted(self.dg_elements): - element = self.dg_elements[facet_dim] - nf = element.space_dimension() - num_facets = len(self.ref_el.get_topology()[facet_dim]) - - # Loop over the number of facets until we find a facet - # with matching dimension and id - for i in range(num_facets): - # Found it! - if facet_dim == entity_dim and i == entity_id: - nonzerovals, = itervalues(element.tabulate(0, points)) - phivals[evalkey][offset:offset+nf, :] = nonzerovals - - offset += nf - - # If asking for gradient evaluations, insert TraceError in - # gradient slots - if order > 0: - msg = "Gradients on trace elements are not well-defined." - for key in phivals: - if key != evalkey: - phivals[key] = TraceError(msg) + # Insert non-zero values in appropriate place + phivals[evalkey][sindex:eindex, :] = nonzerovals return phivals From e061ec06cbe0e2b680f1d2445d57514e23284919 Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Thu, 6 Apr 2017 16:22:06 +0100 Subject: [PATCH 12/12] Use slice notation --- FIAT/hdiv_trace.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/FIAT/hdiv_trace.py b/FIAT/hdiv_trace.py index 8e49fea75..4e400193b 100644 --- a/FIAT/hdiv_trace.py +++ b/FIAT/hdiv_trace.py @@ -217,8 +217,7 @@ def tabulate(self, order, points, entity=None): element = self.dg_elements[facet_sd] nf = element.space_dimension() nonzerovals, = itervalues(element.tabulate(order, new_points)) - sindex = nf * unique_facet - eindex = nf * (unique_facet + 1) + indices = slice(nf * unique_facet, nf * (unique_facet + 1)) else: entity_dim, _ = entity @@ -246,8 +245,7 @@ def tabulate(self, order, points, entity=None): # Found it! Grab insertion indices if (facet_dim, i) == entity: nonzerovals, = itervalues(element.tabulate(0, points)) - sindex = offset - eindex = offset + nf + indices = slice(offset, offset + nf) offset += nf @@ -260,7 +258,7 @@ def tabulate(self, order, points, entity=None): phivals[key] = TraceError(msg) # Insert non-zero values in appropriate place - phivals[evalkey][sindex:eindex, :] = nonzerovals + phivals[evalkey][indices, :] = nonzerovals return phivals