From 7ffda5895e628781a4f0c9bf140d7d2d3dd92812 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Mon, 3 Feb 2020 13:49:34 -0600 Subject: [PATCH 01/10] physical coords --- tsfc/fem.py | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/tsfc/fem.py b/tsfc/fem.py index 5c226384..c76cd623 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -23,6 +23,7 @@ from FIAT.reference_element import make_affine_mapping +from FIAT import ufc_simplex import gem from gem.node import traversal @@ -151,9 +152,19 @@ def jacobian_at(self, point): return map_expr_dag(context.translator, expr) def reference_normals(self): - return gem.Literal(numpy.asarray([self.interface.fiat_cell.compute_normal(i) for i in range(3)])) + if isinstance(self.interface.fiat_cell, FIAT.ufc_simplex) \ + and self.interface.fiat_cell.get_spatial_dimension() == 2: + return gem.Literal(numpy.asarray([self.interface.fiat_cell.compute_normal(i) for i in range(3)])) + else: + raise NotImplementedError("Only works for triangles for now") def physical_tangents(self): + if isinstance(self.interface.fiat_cell, FIAT.ufc_simplex) \ + and self.interface.fiat_cell.get_spatial_dimension() == 2: + return gem.Literal(numpy.asarray([self.interface.fiat_cell.compute_normal(i) for i in range(3)])) + else: + raise NotImplementedError("Only works for triangles for now") + rts = [self.interface.fiat_cell.compute_tangents(1, f)[0] for f in range(3)] jac = self.jacobian_at([1/3, 1/3]) @@ -181,6 +192,24 @@ def physical_edge_lengths(self): context = PointSetContext(**config) return map_expr_dag(context.translator, expr) + def physical_points(self, point_set): + """Converts point_set from reference to physical space""" + expr = SpatialCoordinate(self.mt.terminal.ufl_domain()) + if self.mt.restriction == '+': + expr = PositiveRestricted(expr) + elif self.mt.restriction == '-': + expr = NegativeRestricted(expr) + expr = preprocess_expression(expr) + config = {"point_set": point_set} + config.update(self.config) + context = PointSetContext(**config) + return map_expr_dag(context.translator, expr) + + def physical_vertices(self): + vs = PointSet(self.interface.fiat_cell.vertices) + return self.physical_points(vs) + + def needs_coordinate_mapping(element): """Does this UFL element require a CoordinateMapping for translation?""" From e2ac2748bd8669e27c65c6ca5cadc938ad536623 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Wed, 5 Feb 2020 16:55:52 -0600 Subject: [PATCH 02/10] work on interface --- tsfc/fem.py | 11 ++++++++--- tsfc/finatinterface.py | 1 + 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/tsfc/fem.py b/tsfc/fem.py index c76cd623..959494cf 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -31,7 +31,7 @@ from gem.unconcatenate import unconcatenate from gem.utils import cached_property -from finat.physically_mapped import PhysicalGeometry, PhysicallyMappedElement +from finat.physically_mapped import PhysicalGeometry, NeedsCoordinateMappingElement from finat.point_set import PointSet, PointSingleton from finat.quadrature import make_quadrature @@ -207,7 +207,12 @@ def physical_points(self, point_set): def physical_vertices(self): vs = PointSet(self.interface.fiat_cell.vertices) - return self.physical_points(vs) + pvs = self.physical_points(vs) + return pvs + # print(pvs) + # i = gem.Index() + # j = pvs.free_indices + # return gem.ComponentTensor(gem.Indexed(pvs, (i,)), j + (i,)) @@ -216,7 +221,7 @@ def needs_coordinate_mapping(element): if element.family() == 'Real': return False else: - return isinstance(create_element(element), PhysicallyMappedElement) + return isinstance(create_element(element), NeedsCoordinateMappingElement) class PointSetContext(ContextBase): diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index 84c9fddf..a8cd8457 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -71,6 +71,7 @@ "Discontinuous Lagrange L2": finat.DiscontinuousLagrange, "Gauss-Legendre L2": finat.GaussLegendre, "DQ L2": None, + "Sphys": finat.DirectSerendipity, } """A :class:`.dict` mapping UFL element family names to their FInAT-equivalent constructors. If the value is ``None``, the UFL From 08391b48b3b58f4e437fc42455fff6597f11ab40 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 7 Feb 2020 12:56:42 -0600 Subject: [PATCH 03/10] Work on mapping point sets --- tsfc/fem.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/tsfc/fem.py b/tsfc/fem.py index 959494cf..e3ee0ff0 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -203,16 +203,13 @@ def physical_points(self, point_set): config = {"point_set": point_set} config.update(self.config) context = PointSetContext(**config) - return map_expr_dag(context.translator, expr) - + mapped = map_expr_dag(context.translator, expr) + indices = tuple(gem.Index() for _ in mapped.shape) + return gem.ComponentTensor(gem.Indexed(mapped, indices), point_set.indices + indices) + def physical_vertices(self): vs = PointSet(self.interface.fiat_cell.vertices) - pvs = self.physical_points(vs) - return pvs - # print(pvs) - # i = gem.Index() - # j = pvs.free_indices - # return gem.ComponentTensor(gem.Indexed(pvs, (i,)), j + (i,)) + return self.physical_points(vs) From 02815090326edabf410e79a9a83139d54a1f8e9a Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Sat, 8 Feb 2020 13:03:29 +0000 Subject: [PATCH 04/10] Fix incorrect logic --- tsfc/fem.py | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/tsfc/fem.py b/tsfc/fem.py index e3ee0ff0..3f9ac24d 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -152,17 +152,14 @@ def jacobian_at(self, point): return map_expr_dag(context.translator, expr) def reference_normals(self): - if isinstance(self.interface.fiat_cell, FIAT.ufc_simplex) \ - and self.interface.fiat_cell.get_spatial_dimension() == 2: - return gem.Literal(numpy.asarray([self.interface.fiat_cell.compute_normal(i) for i in range(3)])) - else: + if not (isinstance(self.interface.fiat_cell, ufc_simplex) and + self.interface.fiat_cell.get_spatial_dimension() == 2): raise NotImplementedError("Only works for triangles for now") + return gem.Literal(numpy.asarray([self.interface.fiat_cell.compute_normal(i) for i in range(3)])) def physical_tangents(self): - if isinstance(self.interface.fiat_cell, FIAT.ufc_simplex) \ - and self.interface.fiat_cell.get_spatial_dimension() == 2: - return gem.Literal(numpy.asarray([self.interface.fiat_cell.compute_normal(i) for i in range(3)])) - else: + if not (isinstance(self.interface.fiat_cell, ufc_simplex) and + self.interface.fiat_cell.get_spatial_dimension() == 2): raise NotImplementedError("Only works for triangles for now") rts = [self.interface.fiat_cell.compute_tangents(1, f)[0] for f in range(3)] From e97273536c627272150d243ce6df3b72840429c8 Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Sat, 8 Feb 2020 13:04:15 +0000 Subject: [PATCH 05/10] Add error checking in jacobian_at --- tsfc/fem.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tsfc/fem.py b/tsfc/fem.py index 3f9ac24d..bb6d17e6 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -139,7 +139,9 @@ def cell_size(self): return self.interface.cell_size(self.mt.restriction) def jacobian_at(self, point): + ps = PointSingleton(point) expr = Jacobian(self.mt.terminal.ufl_domain()) + assert ps.expression.shape == (expr.ufl_domain().topological_dimension(), ) if self.mt.restriction == '+': expr = PositiveRestricted(expr) elif self.mt.restriction == '-': From 024220513527d00c7d2ba9f419514363acad4c76 Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Sat, 8 Feb 2020 13:04:26 +0000 Subject: [PATCH 06/10] Implement push-forward of subentity point sets Just need to send the integration_dim and entity_ids through the PointSetContext. --- tsfc/fem.py | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/tsfc/fem.py b/tsfc/fem.py index bb6d17e6..187330d2 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -147,8 +147,7 @@ def jacobian_at(self, point): elif self.mt.restriction == '-': expr = NegativeRestricted(expr) expr = preprocess_expression(expr) - - config = {"point_set": PointSingleton(point)} + config = {"point_set": ps} config.update(self.config) context = PointSetContext(**config) return map_expr_dag(context.translator, expr) @@ -191,9 +190,15 @@ def physical_edge_lengths(self): context = PointSetContext(**config) return map_expr_dag(context.translator, expr) - def physical_points(self, point_set): + def physical_points(self, point_set, entity=None): """Converts point_set from reference to physical space""" expr = SpatialCoordinate(self.mt.terminal.ufl_domain()) + point_shape, = point_set.expression.shape + if entity is not None: + e, _ = entity + assert point_shape == e + else: + assert point_shape == expr.ufl_domain().topological_dimension() if self.mt.restriction == '+': expr = PositiveRestricted(expr) elif self.mt.restriction == '-': @@ -201,17 +206,19 @@ def physical_points(self, point_set): expr = preprocess_expression(expr) config = {"point_set": point_set} config.update(self.config) + if entity is not None: + config.update({name: getattr(self.interface, name) + for name in ["integration_dim", "entity_ids"]}) context = PointSetContext(**config) mapped = map_expr_dag(context.translator, expr) indices = tuple(gem.Index() for _ in mapped.shape) - return gem.ComponentTensor(gem.Indexed(mapped, indices), point_set.indices + indices) - + return gem.ComponentTensor(gem.Indexed(mapped, indices), point_set.indices + indices) + def physical_vertices(self): vs = PointSet(self.interface.fiat_cell.vertices) return self.physical_points(vs) - def needs_coordinate_mapping(element): """Does this UFL element require a CoordinateMapping for translation?""" if element.family() == 'Real': From a1abff57bfb17c775375178e882e8678ad1c2944 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 26 Jun 2020 12:24:06 -0500 Subject: [PATCH 07/10] Rename element --- tsfc/finatinterface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index a8cd8457..96209f59 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -71,7 +71,7 @@ "Discontinuous Lagrange L2": finat.DiscontinuousLagrange, "Gauss-Legendre L2": finat.GaussLegendre, "DQ L2": None, - "Sphys": finat.DirectSerendipity, + "Sdirect": finat.DirectSerendipity, } """A :class:`.dict` mapping UFL element family names to their FInAT-equivalent constructors. If the value is ``None``, the UFL From fd64a007c5ce3783746e3519f4b1302119369e1b Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 3 Jul 2020 09:52:03 -0500 Subject: [PATCH 08/10] Clean up merge rubbish --- tsfc/fem.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/tsfc/fem.py b/tsfc/fem.py index 4a21d483..8a8c90ec 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -147,12 +147,7 @@ def jacobian_at(self, point): expr = PositiveRestricted(expr) elif self.mt.restriction == '-': expr = NegativeRestricted(expr) -# <<<<<<< HEAD -# expr = preprocess_expression(expr) -# config = {"point_set": ps} -# ======= config = {"point_set": PointSingleton(point)} -# >>>>>>> origin/master config.update(self.config) context = PointSetContext(**config) expr = preprocess_expression(expr, complex_mode=context.complex_mode) From b124ed1959bbe6f2f5181fd59190dab7b89bdcf3 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 10 Jul 2020 14:21:14 -0500 Subject: [PATCH 09/10] Fix a type check in CoordinateMapping causing a few tests to fail. --- tsfc/fem.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tsfc/fem.py b/tsfc/fem.py index 8a8c90ec..b9b856f8 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -23,7 +23,7 @@ from FIAT.reference_element import make_affine_mapping -from FIAT import ufc_simplex +from FIAT.reference_element import UFCSimplex import gem from gem.node import traversal @@ -154,13 +154,13 @@ def jacobian_at(self, point): return map_expr_dag(context.translator, expr) def reference_normals(self): - if not (isinstance(self.interface.fiat_cell, ufc_simplex) and + if not (isinstance(self.interface.fiat_cell, UFCSimplex) and self.interface.fiat_cell.get_spatial_dimension() == 2): raise NotImplementedError("Only works for triangles for now") return gem.Literal(numpy.asarray([self.interface.fiat_cell.compute_normal(i) for i in range(3)])) def physical_tangents(self): - if not (isinstance(self.interface.fiat_cell, ufc_simplex) and + if not (isinstance(self.interface.fiat_cell, UFCSimplex) and self.interface.fiat_cell.get_spatial_dimension() == 2): raise NotImplementedError("Only works for triangles for now") From 680293bf5aad6fc53a5fde2236a923740845b571 Mon Sep 17 00:00:00 2001 From: Rob Kirby Date: Fri, 21 Aug 2020 12:25:12 -0500 Subject: [PATCH 10/10] keep up with UFL renaming --- tsfc/finatinterface.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py index 96209f59..2a686eab 100644 --- a/tsfc/finatinterface.py +++ b/tsfc/finatinterface.py @@ -71,7 +71,7 @@ "Discontinuous Lagrange L2": finat.DiscontinuousLagrange, "Gauss-Legendre L2": finat.GaussLegendre, "DQ L2": None, - "Sdirect": finat.DirectSerendipity, + "Direct Serendipity": finat.DirectSerendipity, } """A :class:`.dict` mapping UFL element family names to their FInAT-equivalent constructors. If the value is ``None``, the UFL