From aff37ec047b05d0bc63724da562e2e1fbc88822c Mon Sep 17 00:00:00 2001 From: Thomas Gibson Date: Tue, 13 Sep 2016 18:30:16 +0100 Subject: [PATCH 1/2] add 'entity' argument to FiniteElement.tabulate Aggregates changes from several commits. --- FIAT/bubble.py | 4 +-- FIAT/discontinuous.py | 4 +-- FIAT/enriched.py | 7 +++-- FIAT/finite_element.py | 19 +++++++++++-- FIAT/hdivcurl.py | 8 +++--- FIAT/restricted.py | 4 +-- FIAT/tensor_product.py | 64 +++++++++++++++++++++++++++++++----------- 7 files changed, 78 insertions(+), 32 deletions(-) diff --git a/FIAT/bubble.py b/FIAT/bubble.py index 991328b6a..2c94e5363 100644 --- a/FIAT/bubble.py +++ b/FIAT/bubble.py @@ -73,9 +73,9 @@ def num_sub_elements(self): def space_dimension(self): return self.fsdim - def tabulate(self, order, points): + def tabulate(self, order, points, entity=None): return dict((k, v[self.first_node_index:, :]) - for k, v in self._element.tabulate(order, points).items()) + for k, v in self._element.tabulate(order, points, entity).items()) def value_shape(self): return self._element.value_shape() diff --git a/FIAT/discontinuous.py b/FIAT/discontinuous.py index 00bf9e961..11eaa93be 100644 --- a/FIAT/discontinuous.py +++ b/FIAT/discontinuous.py @@ -76,10 +76,10 @@ def space_dimension(self): "Return the dimension of the finite element space." return self._element.space_dimension() - def tabulate(self, order, points): + def tabulate(self, order, points, entity=None): """Return tabulated values of derivatives up to given order of basis functions at given points.""" - return self._element.tabulate(order, points) + return self._element.tabulate(order, points, entity) def value_shape(self): "Return the value shape of the finite element functions." diff --git a/FIAT/enriched.py b/FIAT/enriched.py index 4da453a19..4e59d5b82 100644 --- a/FIAT/enriched.py +++ b/FIAT/enriched.py @@ -1,5 +1,6 @@ # Copyright (C) 2008 Robert C. Kirby (Texas Tech University) # Copyright (C) 2013 Andrew T. T. McRae +# Modified by Thomas H. Gibson, 2016 # # This file is part of FIAT. # @@ -99,7 +100,7 @@ def space_dimension(self): # number of dofs just adds return self.A.space_dimension() + self.B.space_dimension() - def tabulate(self, order, points): + def tabulate(self, order, points, entity=None): """Return tabulated values of derivatives up to given order of basis functions at given points.""" @@ -109,8 +110,8 @@ def tabulate(self, order, points): Asd = self.A.space_dimension() Bsd = self.B.space_dimension() - Atab = self.A.tabulate(order, points) - Btab = self.B.tabulate(order, points) + Atab = self.A.tabulate(order, points, entity) + Btab = self.B.tabulate(order, points, entity) npoints = len(points) vs = self.A.value_shape() rank = len(vs) # scalar: 0, vector: 1 diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index 38f6f9aa7..a19fe8285 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -17,6 +17,7 @@ # along with FIAT. If not, see . # # Modified by David A. Ham (david.ham@imperial.ac.uk), 2014 +# Modified by Thomas H. Gibson (t.gibson15@imperial.ac.uk), 2016 from __future__ import absolute_import, print_function, division @@ -138,10 +139,22 @@ def space_dimension(self): "Return the dimension of the finite element space." return self.poly_set.get_num_members() - def tabulate(self, order, points): + def tabulate(self, order, points, entity=None): """Return tabulated values of derivatives up to given order of - basis functions at given points.""" - return self.poly_set.tabulate(points, order) + basis functions at given points. + + :arg order: The maximum order of derivative. + :arg points: An iterable of points. + :arg entity: Optional entity indicating which topological entity + of the reference element to tabulate on. If "None", default + cell-wise tabulation is performed. + """ + 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) + return self.poly_set.tabulate(list(map(transform, points)), order) def value_shape(self): "Return the value shape of the finite element functions." diff --git a/FIAT/hdivcurl.py b/FIAT/hdivcurl.py index c8568912c..49f6823a5 100644 --- a/FIAT/hdivcurl.py +++ b/FIAT/hdivcurl.py @@ -53,12 +53,12 @@ def value_shape(self): # redefine tabulate newelement.old_tabulate = newelement.tabulate - def tabulate(self, order, points): + def tabulate(self, order, points, entity=None): """Return tabulated values of derivatives up to given order of basis functions at given points.""" # don't duplicate what the old function does fine... - old_result = self.old_tabulate(order, points) + old_result = self.old_tabulate(order, points, entity) new_result = {} sd = self.get_reference_element().get_spatial_dimension() for alpha in old_result.keys(): @@ -175,12 +175,12 @@ def value_shape(self): # redefine tabulate newelement.old_tabulate = newelement.tabulate - def tabulate(self, order, points): + def tabulate(self, order, points, entity=None): """Return tabulated values of derivatives up to given order of basis functions at given points.""" # don't duplicate what the old function does fine... - old_result = self.old_tabulate(order, points) + old_result = self.old_tabulate(order, points, entity) new_result = {} sd = self.get_reference_element().get_spatial_dimension() for alpha in old_result.keys(): diff --git a/FIAT/restricted.py b/FIAT/restricted.py index 6f01c300d..45f21552d 100644 --- a/FIAT/restricted.py +++ b/FIAT/restricted.py @@ -85,8 +85,8 @@ def mapping(self): mappings = self._element.mapping() return [mappings[i] for i in self._indices] - def tabulate(self, order, points): - result = self._element.tabulate(order, points) + def tabulate(self, order, points, entity=None): + result = self._element.tabulate(order, points, entity) extracted = {} for (dtuple, values) in sorted_by_key(result): extracted[dtuple] = numpy.array([values[i] for i in self._indices]) diff --git a/FIAT/tensor_product.py b/FIAT/tensor_product.py index fac3414b7..6b7698234 100644 --- a/FIAT/tensor_product.py +++ b/FIAT/tensor_product.py @@ -1,5 +1,6 @@ # Copyright (C) 2008 Robert C. Kirby (Texas Tech University) # Copyright (C) 2013 Andrew T. T. McRae +# Modified by Thomas H. Gibson, 2016 # # This file is part of FIAT. # @@ -240,16 +241,47 @@ def space_dimension(self): # number of dofs just multiplies return self.A.space_dimension() * self.B.space_dimension() - def tabulate(self, order, points): + def tabulate(self, order, points, entity=None): """Return tabulated values of derivatives up to given order of basis functions at given points.""" + Adim = self.A.ref_el.get_spatial_dimension() + Bdim = self.B.ref_el.get_spatial_dimension() + + if entity is None: + entity = (self.ref_el.get_dimension(), 0) + + shape = tuple(len(c.get_topology()[d]) + for c, d in zip(self.ref_el.cells, entity[0])) + idA, idB = numpy.unravel_index(entity[1], shape) + + # Factor the entity argument to get entities of the component elements + entityA_dim, entityB_dim = entity[0] + entityA = (entityA_dim, idA) + entityB = (entityB_dim, idB) + + pointsAdim, pointsBdim = [c.get_spatial_dimension() for c in self.ref_el.construct_subelement(entity[0]).cells] + # Sum dimensions in case the sub-elements are themselves tensor products. + Aslice = slice(Adim) + Bslice = slice(Adim, Adim + Bdim) + pointsA = [point[:pointsAdim] for point in points] + pointsB = [point[pointsAdim:pointsAdim + pointsBdim] for point in points] + # Note that for entities other than cells, the following + # tabulations are already appropriately zero-padded so no + # additional zero padding is required. + try: + Atab = self.A.tabulate(order, pointsA, entityA) + except TraceError as e: + if self.ref_el.get_spatial_dimension() == entityA[0] + entityB[0]: + raise + Atab = e.zeros + + try: + Btab = self.B.tabulate(order, pointsB, entityB) + except TraceError as e: + if self.ref_el.get_spatial_dimension() == entityA[0] + entityB[0]: + raise + Btab = e.zeros - Asdim = self.A.get_reference_element().get_spatial_dimension() - Bsdim = self.B.get_reference_element().get_spatial_dimension() - pointsA = [point[0:Asdim] for point in points] - pointsB = [point[Asdim:Asdim + Bsdim] for point in points] - Atab = self.A.tabulate(order, pointsA) - Btab = self.B.tabulate(order, pointsB) npoints = len(points) # allow 2 scalar-valued FE spaces, or 1 scalar-valued, @@ -266,7 +298,7 @@ def tabulate(self, order, points): raise NotImplementedError("tabulate does not support two vector-valued inputs") result = {} for i in range(order + 1): - alphas = mis(Asdim+Bsdim, i) # thanks, Rob! + alphas = mis(self.ref_el.get_spatial_dimension(), i) # thanks, Rob! for alpha in alphas: if A_valuedim == 0 and B_valuedim == 0: # for each point, get outer product of (A's basis @@ -280,8 +312,8 @@ def tabulate(self, order, points): # Transpose this to get temp[bf][point], # and we are done. temp = numpy.array([numpy.outer( - Atab[alpha[0:Asdim]][..., j], - Btab[alpha[Asdim:Asdim+Bsdim]][..., j]) + Atab[alpha[Aslice]][..., j], + Btab[alpha[Bslice]][..., j]) .ravel() for j in range(npoints)]) result[alpha] = temp.transpose() elif A_valuedim == 1 and B_valuedim == 0: @@ -298,8 +330,8 @@ def tabulate(self, order, points): # finally, transpose the first and last indices # to get temp[bf_i][x/y][point], and we are done. temp = numpy.array([numpy.outer( - Atab[alpha[0:Asdim]][..., j], - Btab[alpha[Asdim:Asdim+Bsdim]][..., j]) + Atab[alpha[Aslice]][..., j], + Btab[alpha[Bslice]][..., j]) for j in range(npoints)]) assert temp.shape[1] % 2 == 0 temp2 = temp.reshape((temp.shape[0], @@ -318,10 +350,10 @@ def tabulate(self, order, points): # flatten middle: temp[point][full bf_i][x/y] # transpose to temp[bf_i][x/y][point] temp = numpy.array([numpy.outer( - Atab[alpha[0:Asdim]][..., j], - Btab[alpha[Asdim:Asdim+Bsdim]][..., j]) - for j in range(len(Atab[alpha[0:Asdim]][0]))]) - assert temp.shape[2] % 2 == 0 + Atab[alpha[Aslice]][..., j], + Btab[alpha[Bslice]][..., j]) + for j in range(len(Atab[alpha[Aslice]][0]))]) + temp2 = temp.reshape((temp.shape[0], temp.shape[1], temp.shape[2]//2, 2))\ .reshape((temp.shape[0], -1, 2))\ From 1e50090ac146a99774ffe84bc26d1b5a59892a3d Mon Sep 17 00:00:00 2001 From: Miklos Homolya Date: Mon, 19 Sep 2016 14:19:56 +0100 Subject: [PATCH 2/2] clean up TensorProductElement.tabulate ... and improve the docstring of FiniteElement.tabulate --- FIAT/finite_element.py | 7 +++--- FIAT/tensor_product.py | 53 ++++++++++++++++-------------------------- 2 files changed, 24 insertions(+), 36 deletions(-) diff --git a/FIAT/finite_element.py b/FIAT/finite_element.py index a19fe8285..bd164ff33 100644 --- a/FIAT/finite_element.py +++ b/FIAT/finite_element.py @@ -145,9 +145,10 @@ def tabulate(self, order, points, entity=None): :arg order: The maximum order of derivative. :arg points: An iterable of points. - :arg entity: Optional entity indicating which topological entity - of the reference element to tabulate on. If "None", default - cell-wise tabulation is performed. + :arg entity: Optional (dimension, entity number) pair + indicating which topological entity of the + reference element to tabulate on. If ``None``, + default cell-wise tabulation is performed. """ if entity is None: entity = (self.ref_el.get_spatial_dimension(), 0) diff --git a/FIAT/tensor_product.py b/FIAT/tensor_product.py index 6b7698234..9de798478 100644 --- a/FIAT/tensor_product.py +++ b/FIAT/tensor_product.py @@ -244,44 +244,31 @@ def space_dimension(self): def tabulate(self, order, points, entity=None): """Return tabulated values of derivatives up to given order of basis functions at given points.""" - Adim = self.A.ref_el.get_spatial_dimension() - Bdim = self.B.ref_el.get_spatial_dimension() - if entity is None: entity = (self.ref_el.get_dimension(), 0) + entity_dim, entity_id = entity shape = tuple(len(c.get_topology()[d]) - for c, d in zip(self.ref_el.cells, entity[0])) - idA, idB = numpy.unravel_index(entity[1], shape) + for c, d in zip(self.ref_el.cells, entity_dim)) + idA, idB = numpy.unravel_index(entity_id, shape) # Factor the entity argument to get entities of the component elements - entityA_dim, entityB_dim = entity[0] + entityA_dim, entityB_dim = entity_dim entityA = (entityA_dim, idA) entityB = (entityB_dim, idB) - pointsAdim, pointsBdim = [c.get_spatial_dimension() for c in self.ref_el.construct_subelement(entity[0]).cells] - # Sum dimensions in case the sub-elements are themselves tensor products. - Aslice = slice(Adim) - Bslice = slice(Adim, Adim + Bdim) + pointsAdim, pointsBdim = [c.get_spatial_dimension() + for c in self.ref_el.construct_subelement(entity_dim).cells] pointsA = [point[:pointsAdim] for point in points] pointsB = [point[pointsAdim:pointsAdim + pointsBdim] for point in points] + + Asdim = self.A.ref_el.get_spatial_dimension() + Bsdim = self.B.ref_el.get_spatial_dimension() # Note that for entities other than cells, the following # tabulations are already appropriately zero-padded so no # additional zero padding is required. - try: - Atab = self.A.tabulate(order, pointsA, entityA) - except TraceError as e: - if self.ref_el.get_spatial_dimension() == entityA[0] + entityB[0]: - raise - Atab = e.zeros - - try: - Btab = self.B.tabulate(order, pointsB, entityB) - except TraceError as e: - if self.ref_el.get_spatial_dimension() == entityA[0] + entityB[0]: - raise - Btab = e.zeros - + Atab = self.A.tabulate(order, pointsA, entityA) + Btab = self.B.tabulate(order, pointsB, entityB) npoints = len(points) # allow 2 scalar-valued FE spaces, or 1 scalar-valued, @@ -298,7 +285,7 @@ def tabulate(self, order, points, entity=None): raise NotImplementedError("tabulate does not support two vector-valued inputs") result = {} for i in range(order + 1): - alphas = mis(self.ref_el.get_spatial_dimension(), i) # thanks, Rob! + alphas = mis(Asdim+Bsdim, i) # thanks, Rob! for alpha in alphas: if A_valuedim == 0 and B_valuedim == 0: # for each point, get outer product of (A's basis @@ -312,8 +299,8 @@ def tabulate(self, order, points, entity=None): # Transpose this to get temp[bf][point], # and we are done. temp = numpy.array([numpy.outer( - Atab[alpha[Aslice]][..., j], - Btab[alpha[Bslice]][..., j]) + Atab[alpha[0:Asdim]][..., j], + Btab[alpha[Asdim:Asdim+Bsdim]][..., j]) .ravel() for j in range(npoints)]) result[alpha] = temp.transpose() elif A_valuedim == 1 and B_valuedim == 0: @@ -330,8 +317,8 @@ def tabulate(self, order, points, entity=None): # finally, transpose the first and last indices # to get temp[bf_i][x/y][point], and we are done. temp = numpy.array([numpy.outer( - Atab[alpha[Aslice]][..., j], - Btab[alpha[Bslice]][..., j]) + Atab[alpha[0:Asdim]][..., j], + Btab[alpha[Asdim:Asdim+Bsdim]][..., j]) for j in range(npoints)]) assert temp.shape[1] % 2 == 0 temp2 = temp.reshape((temp.shape[0], @@ -350,10 +337,10 @@ def tabulate(self, order, points, entity=None): # flatten middle: temp[point][full bf_i][x/y] # transpose to temp[bf_i][x/y][point] temp = numpy.array([numpy.outer( - Atab[alpha[Aslice]][..., j], - Btab[alpha[Bslice]][..., j]) - for j in range(len(Atab[alpha[Aslice]][0]))]) - + Atab[alpha[0:Asdim]][..., j], + Btab[alpha[Asdim:Asdim+Bsdim]][..., j]) + for j in range(len(Atab[alpha[0:Asdim]][0]))]) + assert temp.shape[2] % 2 == 0 temp2 = temp.reshape((temp.shape[0], temp.shape[1], temp.shape[2]//2, 2))\ .reshape((temp.shape[0], -1, 2))\