diff --git a/tsfc/fem.py b/tsfc/fem.py index f70df954..b9b856f8 100644 --- a/tsfc/fem.py +++ b/tsfc/fem.py @@ -23,6 +23,7 @@ from FIAT.reference_element import make_affine_mapping +from FIAT.reference_element import UFCSimplex import gem from gem.node import traversal @@ -30,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 @@ -139,7 +140,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 == '-': @@ -151,9 +154,16 @@ def jacobian_at(self, point): return map_expr_dag(context.translator, expr) def reference_normals(self): + 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, UFCSimplex) 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)] jac = self.jacobian_at([1/3, 1/3]) @@ -181,13 +191,41 @@ def physical_edge_lengths(self): expr = preprocess_expression(expr, complex_mode=context.complex_mode) return map_expr_dag(context.translator, expr) + 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 == '-': + expr = NegativeRestricted(expr) + 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) + + 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': 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 39734b04..81464585 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, + "Direct Serendipity": finat.DirectSerendipity, } """A :class:`.dict` mapping UFL element family names to their FInAT-equivalent constructors. If the value is ``None``, the UFL