From b5f77b5d720544048beadff871ddae24aa36993f Mon Sep 17 00:00:00 2001 From: Jack Betteridge Date: Wed, 23 Oct 2024 15:22:47 +0100 Subject: [PATCH] Remove TSFC --- setup.py | 2 +- tsfc/__init__.py | 22 - tsfc/coffee_mode.py | 81 --- tsfc/driver.py | 334 ----------- tsfc/fem.py | 729 ----------------------- tsfc/finatinterface.py | 359 ----------- tsfc/kernel_args.py | 54 -- tsfc/kernel_interface/__init__.py | 47 -- tsfc/kernel_interface/common.py | 540 ----------------- tsfc/kernel_interface/firedrake_loopy.py | 435 -------------- tsfc/logging.py | 6 - tsfc/loopy.py | 549 ----------------- tsfc/modified_terminals.py | 179 ------ tsfc/parameters.py | 36 -- tsfc/spectral.py | 193 ------ tsfc/tensor.py | 93 --- tsfc/ufl2gem.py | 160 ----- tsfc/ufl_utils.py | 485 --------------- tsfc/vanilla.py | 47 -- 19 files changed, 1 insertion(+), 4350 deletions(-) delete mode 100644 tsfc/__init__.py delete mode 100644 tsfc/coffee_mode.py delete mode 100644 tsfc/driver.py delete mode 100644 tsfc/fem.py delete mode 100644 tsfc/finatinterface.py delete mode 100644 tsfc/kernel_args.py delete mode 100644 tsfc/kernel_interface/__init__.py delete mode 100644 tsfc/kernel_interface/common.py delete mode 100644 tsfc/kernel_interface/firedrake_loopy.py delete mode 100644 tsfc/logging.py delete mode 100644 tsfc/loopy.py delete mode 100644 tsfc/modified_terminals.py delete mode 100644 tsfc/parameters.py delete mode 100644 tsfc/spectral.py delete mode 100644 tsfc/tensor.py delete mode 100644 tsfc/ufl2gem.py delete mode 100644 tsfc/ufl_utils.py delete mode 100644 tsfc/vanilla.py diff --git a/setup.py b/setup.py index e300a605..bb238492 100644 --- a/setup.py +++ b/setup.py @@ -4,4 +4,4 @@ setup(name="tsfc", version=version, - packages=["gem", "tsfc", "tsfc.kernel_interface"]) + packages=["gem"]) diff --git a/tsfc/__init__.py b/tsfc/__init__.py deleted file mode 100644 index f9075c71..00000000 --- a/tsfc/__init__.py +++ /dev/null @@ -1,22 +0,0 @@ -from tsfc.driver import compile_form, compile_expression_dual_evaluation # noqa: F401 -from tsfc.parameters import default_parameters # noqa: F401 - -try: - from firedrake_citations import Citations - Citations().add("Kirby2006", """ -@Article{Kirby2006, - author = {Kirby, Robert C. and Logg, Anders}, - title = {A Compiler for Variational Forms}, - journal = {ACM Trans. Math. Softw.}, - year = 2006, - volume = 32, - number = 3, - pages = {417--444}, - month = sep, - numpages = 28, - doi = {10.1145/1163641.1163644}, - acmid = 1163644, -}""") - del Citations -except ImportError: - pass diff --git a/tsfc/coffee_mode.py b/tsfc/coffee_mode.py deleted file mode 100644 index 632b915b..00000000 --- a/tsfc/coffee_mode.py +++ /dev/null @@ -1,81 +0,0 @@ -from functools import partial, reduce - -from gem.node import traversal, Memoizer -from gem.gem import Failure, Sum, index_sum -from gem.optimise import replace_division, unroll_indexsum -from gem.refactorise import collect_monomials -from gem.unconcatenate import unconcatenate -from gem.coffee import optimise_monomial_sum -from gem.utils import groupby - -import tsfc.spectral as spectral - - -def Integrals(expressions, quadrature_multiindex, argument_multiindices, parameters): - """Constructs an integral representation for each GEM integrand - expression. - - :arg expressions: integrand multiplied with quadrature weight; - multi-root GEM expression DAG - :arg quadrature_multiindex: quadrature multiindex (tuple) - :arg argument_multiindices: tuple of argument multiindices, - one multiindex for each argument - :arg parameters: parameters dictionary - - :returns: list of integral representations - """ - # Unroll - max_extent = parameters["unroll_indexsum"] - if max_extent: - def predicate(index): - return index.extent <= max_extent - expressions = unroll_indexsum(expressions, predicate=predicate) - # Integral representation: just a GEM expression - return replace_division([index_sum(e, quadrature_multiindex) for e in expressions]) - - -def flatten(var_reps, index_cache): - """Flatten mode-specific intermediate representation to a series of - assignments. - - :arg var_reps: series of (return variable, [integral representation]) pairs - :arg index_cache: cache :py:class:`dict` for :py:func:`unconcatenate` - - :returns: series of (return variable, GEM expression root) pairs - """ - assignments = unconcatenate([(variable, reduce(Sum, reps)) - for variable, reps in var_reps], - cache=index_cache) - - def group_key(assignment): - variable, expression = assignment - return variable.free_indices - - for argument_indices, assignment_group in groupby(assignments, group_key): - variables, expressions = zip(*assignment_group) - expressions = optimise_expressions(expressions, argument_indices) - for var, expr in zip(variables, expressions): - yield (var, expr) - - -finalise_options = dict(remove_componenttensors=False) - - -def optimise_expressions(expressions, argument_indices): - """Perform loop optimisations on GEM DAGs - - :arg expressions: list of GEM DAGs - :arg argument_indices: tuple of argument indices - - :returns: list of optimised GEM DAGs - """ - # Skip optimisation for if Failure node is present - for n in traversal(expressions): - if isinstance(n, Failure): - return expressions - - # Apply argument factorisation unconditionally - classifier = partial(spectral.classify, set(argument_indices), - delta_inside=Memoizer(spectral._delta_inside)) - monomial_sums = collect_monomials(expressions, classifier) - return [optimise_monomial_sum(ms, argument_indices) for ms in monomial_sums] diff --git a/tsfc/driver.py b/tsfc/driver.py deleted file mode 100644 index 95fd43fb..00000000 --- a/tsfc/driver.py +++ /dev/null @@ -1,334 +0,0 @@ -import collections -import time -import sys -from itertools import chain -from finat.physically_mapped import DirectlyDefinedElement, PhysicallyMappedElement - -import ufl -from ufl.algorithms import extract_arguments, extract_coefficients -from ufl.algorithms.analysis import has_type -from ufl.classes import Form, GeometricQuantity -from ufl.domain import extract_unique_domain - -import gem -import gem.impero_utils as impero_utils - -import finat - -from tsfc import fem, ufl_utils -from tsfc.logging import logger -from tsfc.parameters import default_parameters, is_complex -from tsfc.ufl_utils import apply_mapping, extract_firedrake_constants -import tsfc.kernel_interface.firedrake_loopy as firedrake_interface_loopy - -# To handle big forms. The various transformations might need a deeper stack -sys.setrecursionlimit(3000) - - -TSFCIntegralDataInfo = collections.namedtuple("TSFCIntegralDataInfo", - ["domain", "integral_type", "subdomain_id", "domain_number", - "arguments", - "coefficients", "coefficient_numbers"]) -TSFCIntegralDataInfo.__doc__ = """ - Minimal set of objects for kernel builders. - - domain - The mesh. - integral_type - The type of integral. - subdomain_id - What is the subdomain id for this kernel. - domain_number - Which domain number in the original form - does this kernel correspond to (can be used to index into - original_form.ufl_domains() to get the correct domain). - coefficients - A list of coefficients. - coefficient_numbers - A list of which coefficients from the - form the kernel needs. - - This is a minimal set of objects that kernel builders need to - construct a kernel from :attr:`integrals` of :class:`~ufl.IntegralData`. - """ - - -def compile_form(form, prefix="form", parameters=None, interface=None, diagonal=False, log=False): - """Compiles a UFL form into a set of assembly kernels. - - :arg form: UFL form - :arg prefix: kernel name will start with this string - :arg parameters: parameters object - :arg diagonal: Are we building a kernel for the diagonal of a rank-2 element tensor? - :arg log: bool if the Kernel should be profiled with Log events - :returns: list of kernels - """ - cpu_time = time.time() - - assert isinstance(form, Form) - - GREEN = "\033[1;37;32m%s\033[0m" - - # Determine whether in complex mode: - complex_mode = parameters and is_complex(parameters.get("scalar_type")) - fd = ufl_utils.compute_form_data(form, complex_mode=complex_mode) - logger.info(GREEN % "compute_form_data finished in %g seconds.", time.time() - cpu_time) - - kernels = [] - for integral_data in fd.integral_data: - start = time.time() - kernel = compile_integral(integral_data, fd, prefix, parameters, interface=interface, diagonal=diagonal, log=log) - if kernel is not None: - kernels.append(kernel) - logger.info(GREEN % "compile_integral finished in %g seconds.", time.time() - start) - - logger.info(GREEN % "TSFC finished in %g seconds.", time.time() - cpu_time) - return kernels - - -def compile_integral(integral_data, form_data, prefix, parameters, interface, *, diagonal=False, log=False): - """Compiles a UFL integral into an assembly kernel. - - :arg integral_data: UFL integral data - :arg form_data: UFL form data - :arg prefix: kernel name will start with this string - :arg parameters: parameters object - :arg interface: backend module for the kernel interface - :arg diagonal: Are we building a kernel for the diagonal of a rank-2 element tensor? - :arg log: bool if the Kernel should be profiled with Log events - :returns: a kernel constructed by the kernel interface - """ - if integral_data.domain.ufl_cell().cellname() == "hexahedron" and \ - integral_data.integral_type == "interior_facet": - raise NotImplementedError("interior facet integration in hex meshes not currently supported") - parameters = preprocess_parameters(parameters) - if interface is None: - interface = firedrake_interface_loopy.KernelBuilder - scalar_type = parameters["scalar_type"] - integral_type = integral_data.integral_type - if integral_type.startswith("interior_facet") and diagonal: - raise NotImplementedError("Sorry, we can't assemble the diagonal of a form for interior facet integrals") - mesh = integral_data.domain - arguments = form_data.preprocessed_form.arguments() - kernel_name = f"{prefix}_{integral_type}_integral" - # Dict mapping domains to index in original_form.ufl_domains() - domain_numbering = form_data.original_form.domain_numbering() - domain_number = domain_numbering[integral_data.domain] - coefficients = [form_data.function_replace_map[c] for c in integral_data.integral_coefficients] - # This is which coefficient in the original form the - # current coefficient is. - # Consider f*v*dx + g*v*ds, the full form contains two - # coefficients, but each integral only requires one. - coefficient_numbers = tuple(form_data.original_coefficient_positions[i] - for i, (_, enabled) in enumerate(zip(form_data.reduced_coefficients, integral_data.enabled_coefficients)) - if enabled) - integral_data_info = TSFCIntegralDataInfo(domain=integral_data.domain, - integral_type=integral_data.integral_type, - subdomain_id=integral_data.subdomain_id, - domain_number=domain_number, - arguments=arguments, - coefficients=coefficients, - coefficient_numbers=coefficient_numbers) - builder = interface(integral_data_info, - scalar_type, - diagonal=diagonal) - builder.set_coordinates(mesh) - builder.set_cell_sizes(mesh) - builder.set_coefficients(integral_data, form_data) - # TODO: We do not want pass constants to kernels that do not need them - # so we should attach the constants to integral data instead - builder.set_constants(form_data.constants) - ctx = builder.create_context() - for integral in integral_data.integrals: - params = parameters.copy() - params.update(integral.metadata()) # integral metadata overrides - integrand = ufl.replace(integral.integrand(), form_data.function_replace_map) - integrand_exprs = builder.compile_integrand(integrand, params, ctx) - integral_exprs = builder.construct_integrals(integrand_exprs, params) - builder.stash_integrals(integral_exprs, params, ctx) - return builder.construct_kernel(kernel_name, ctx, log) - - -def preprocess_parameters(parameters): - if parameters is None: - parameters = default_parameters() - else: - _ = default_parameters() - _.update(parameters) - parameters = _ - # Remove these here, they're handled later on. - if parameters.get("quadrature_degree") in ["auto", "default", None, -1, "-1"]: - del parameters["quadrature_degree"] - if parameters.get("quadrature_rule") in ["auto", "default", None]: - del parameters["quadrature_rule"] - return parameters - - -def compile_expression_dual_evaluation(expression, to_element, ufl_element, *, - domain=None, interface=None, - parameters=None, log=False): - """Compile a UFL expression to be evaluated against a compile-time known reference element's dual basis. - - Useful for interpolating UFL expressions into e.g. N1curl spaces. - - :arg expression: UFL expression - :arg to_element: A FInAT element for the target space - :arg ufl_element: The UFL element of the target space. - :arg domain: optional UFL domain the expression is defined on (required when expression contains no domain). - :arg interface: backend module for the kernel interface - :arg parameters: parameters object - :arg log: bool if the Kernel should be profiled with Log events - :returns: Loopy-based ExpressionKernel object. - """ - if parameters is None: - parameters = default_parameters() - else: - _ = default_parameters() - _.update(parameters) - parameters = _ - - # Determine whether in complex mode - complex_mode = is_complex(parameters["scalar_type"]) - - if isinstance(to_element, (PhysicallyMappedElement, DirectlyDefinedElement)): - raise NotImplementedError("Don't know how to interpolate onto zany spaces, sorry") - - orig_expression = expression - - # Map into reference space - expression = apply_mapping(expression, ufl_element, domain) - - # Apply UFL preprocessing - expression = ufl_utils.preprocess_expression(expression, - complex_mode=complex_mode) - - # Initialise kernel builder - if interface is None: - # Delayed import, loopy is a runtime dependency - from tsfc.kernel_interface.firedrake_loopy import ExpressionKernelBuilder as interface - - builder = interface(parameters["scalar_type"]) - arguments = extract_arguments(expression) - argument_multiindices = tuple(builder.create_element(arg.ufl_element()).get_indices() - for arg in arguments) - - # Replace coordinates (if any) unless otherwise specified by kwarg - if domain is None: - domain = extract_unique_domain(expression) - assert domain is not None - - # Collect required coefficients and determine numbering - coefficients = extract_coefficients(expression) - orig_coefficients = extract_coefficients(orig_expression) - coefficient_numbers = tuple(orig_coefficients.index(c) for c in coefficients) - builder.set_coefficient_numbers(coefficient_numbers) - - needs_external_coords = False - if has_type(expression, GeometricQuantity) or any(fem.needs_coordinate_mapping(c.ufl_element()) for c in coefficients): - # Create a fake coordinate coefficient for a domain. - coords_coefficient = ufl.Coefficient(ufl.FunctionSpace(domain, domain.ufl_coordinate_element())) - builder.domain_coordinate[domain] = coords_coefficient - builder.set_cell_sizes(domain) - coefficients = [coords_coefficient] + coefficients - needs_external_coords = True - builder.set_coefficients(coefficients) - - constants = extract_firedrake_constants(expression) - builder.set_constants(constants) - - # Split mixed coefficients - expression = ufl_utils.split_coefficients(expression, builder.coefficient_split) - - # Set up kernel config for translation of UFL expression to gem - kernel_cfg = dict(interface=builder, - ufl_cell=domain.ufl_cell(), - # FIXME: change if we ever implement - # interpolation on facets. - integral_type="cell", - argument_multiindices=argument_multiindices, - index_cache={}, - scalar_type=parameters["scalar_type"]) - - # Allow interpolation onto QuadratureElements to refer to the quadrature - # rule they represent - if isinstance(to_element, finat.QuadratureElement): - kernel_cfg["quadrature_rule"] = to_element._rule - - # Create callable for translation of UFL expression to gem - fn = DualEvaluationCallable(expression, kernel_cfg) - - # Get the gem expression for dual evaluation and corresponding basis - # indices needed for compilation of the expression - evaluation, basis_indices = to_element.dual_evaluation(fn) - - # Build kernel body - return_indices = basis_indices + tuple(chain(*argument_multiindices)) - return_shape = tuple(i.extent for i in return_indices) - return_var = gem.Variable('A', return_shape) - return_expr = gem.Indexed(return_var, return_indices) - - # TODO: one should apply some GEM optimisations as in assembly, - # but we don't for now. - evaluation, = impero_utils.preprocess_gem([evaluation]) - impero_c = impero_utils.compile_gem([(return_expr, evaluation)], return_indices) - index_names = dict((idx, "p%d" % i) for (i, idx) in enumerate(basis_indices)) - # Handle kernel interface requirements - builder.register_requirements([evaluation]) - builder.set_output(return_var) - # Build kernel tuple - return builder.construct_kernel(impero_c, index_names, needs_external_coords, log=log) - - -class DualEvaluationCallable(object): - """ - Callable representing a function to dual evaluate. - - When called, this takes in a - :class:`finat.point_set.AbstractPointSet` and returns a GEM - expression for evaluation of the function at those points. - - :param expression: UFL expression for the function to dual evaluate. - :param kernel_cfg: A kernel configuration for creation of a - :class:`GemPointContext` or a :class:`PointSetContext` - - Not intended for use outside of - :func:`compile_expression_dual_evaluation`. - """ - def __init__(self, expression, kernel_cfg): - self.expression = expression - self.kernel_cfg = kernel_cfg - - def __call__(self, ps): - """The function to dual evaluate. - - :param ps: The :class:`finat.point_set.AbstractPointSet` for - evaluating at - :returns: a gem expression representing the evaluation of the - input UFL expression at the given point set ``ps``. - For point set points with some shape ``(*value_shape)`` - (i.e. ``()`` for scalar points ``(x)`` for vector points - ``(x, y)`` for tensor points etc) then the gem expression - has shape ``(*value_shape)`` and free indices corresponding - to the input :class:`finat.point_set.AbstractPointSet`'s - free indices alongside any input UFL expression free - indices. - """ - - if not isinstance(ps, finat.point_set.AbstractPointSet): - raise ValueError("Callable argument not a point set!") - - # Avoid modifying saved kernel config - kernel_cfg = self.kernel_cfg.copy() - - if isinstance(ps, finat.point_set.UnknownPointSet): - # Run time known points - kernel_cfg.update(point_indices=ps.indices, point_expr=ps.expression) - # GemPointContext's aren't allowed to have quadrature rules - kernel_cfg.pop("quadrature_rule", None) - translation_context = fem.GemPointContext(**kernel_cfg) - else: - # Compile time known points - kernel_cfg.update(point_set=ps) - translation_context = fem.PointSetContext(**kernel_cfg) - - gem_expr, = fem.compile_ufl(self.expression, translation_context, point_sum=False) - # In some cases ps.indices may be dropped from expr, but nothing - # new should now appear - argument_multiindices = kernel_cfg["argument_multiindices"] - assert set(gem_expr.free_indices) <= set(chain(ps.indices, *argument_multiindices)) - - return gem_expr diff --git a/tsfc/fem.py b/tsfc/fem.py deleted file mode 100644 index 27439f1c..00000000 --- a/tsfc/fem.py +++ /dev/null @@ -1,729 +0,0 @@ -"""Functions to translate UFL finite element objects and reference -geometric quantities into GEM expressions.""" - -import collections -import itertools -from functools import singledispatch - -import gem -import numpy -import ufl -from FIAT.reference_element import UFCSimplex, make_affine_mapping -from finat.physically_mapped import (NeedsCoordinateMappingElement, - PhysicalGeometry) -from finat.point_set import PointSet, PointSingleton -from finat.quadrature import make_quadrature -from gem.node import traversal -from gem.optimise import constant_fold_zero, ffc_rounding -from gem.unconcatenate import unconcatenate -from gem.utils import cached_property -from ufl.classes import (Argument, CellCoordinate, CellEdgeVectors, - CellFacetJacobian, CellOrientation, CellOrigin, - CellVertices, CellVolume, Coefficient, FacetArea, - FacetCoordinate, GeometricQuantity, Jacobian, - JacobianDeterminant, NegativeRestricted, - PositiveRestricted, QuadratureWeight, - ReferenceCellEdgeVectors, ReferenceCellVolume, - ReferenceFacetVolume, ReferenceNormal, - SpatialCoordinate) -from ufl.corealg.map_dag import map_expr_dag, map_expr_dags -from ufl.corealg.multifunction import MultiFunction -from ufl.domain import extract_unique_domain - -from tsfc import ufl2gem -from tsfc.finatinterface import as_fiat_cell, create_element -from tsfc.kernel_interface import ProxyKernelInterface -from tsfc.modified_terminals import (analyse_modified_terminal, - construct_modified_terminal) -from tsfc.parameters import is_complex -from tsfc.ufl_utils import (ModifiedTerminalMixin, PickRestriction, - TSFCConstantMixin, entity_avg, one_times, - preprocess_expression, simplify_abs) - - -class ContextBase(ProxyKernelInterface): - """Common UFL -> GEM translation context.""" - - keywords = ('ufl_cell', - 'fiat_cell', - 'integral_type', - 'integration_dim', - 'entity_ids', - 'argument_multiindices', - 'facetarea', - 'index_cache', - 'scalar_type') - - def __init__(self, interface, **kwargs): - ProxyKernelInterface.__init__(self, interface) - - invalid_keywords = set(kwargs.keys()) - set(self.keywords) - if invalid_keywords: - raise ValueError("unexpected keyword argument '{0}'".format(invalid_keywords.pop())) - self.__dict__.update(kwargs) - - @cached_property - def fiat_cell(self): - return as_fiat_cell(self.ufl_cell) - - @cached_property - def integration_dim(self): - return self.fiat_cell.get_dimension() - - entity_ids = [0] - - @cached_property - def epsilon(self): - return numpy.finfo(self.scalar_type).resolution - - @cached_property - def complex_mode(self): - return is_complex(self.scalar_type) - - def entity_selector(self, callback, restriction): - """Selects code for the correct entity at run-time. Callback - generates code for a specified entity. - - This function passes ``callback`` the entity number. - - :arg callback: A function to be called with an entity number - that generates code for that entity. - :arg restriction: Restriction of the modified terminal, used - for entity selection. - """ - if len(self.entity_ids) == 1: - return callback(self.entity_ids[0]) - else: - f = self.entity_number(restriction) - return gem.select_expression(list(map(callback, self.entity_ids)), f) - - argument_multiindices = () - - @cached_property - def index_cache(self): - return {} - - @cached_property - def translator(self): - # NOTE: reference cycle! - return Translator(self) - - -class CoordinateMapping(PhysicalGeometry): - """Callback class that provides physical geometry to FInAT elements. - - Required for elements whose basis transformation requires physical - geometry such as Argyris and Hermite. - - :arg mt: The modified terminal whose element will be tabulated. - :arg interface: The interface carrying information (generally a - :class:`PointSetContext`). - """ - def __init__(self, mt, interface): - super().__init__() - self.mt = mt - self.interface = interface - - def preprocess(self, expr, context): - """Preprocess a UFL expression for translation. - - :arg expr: A UFL expression - :arg context: The translation context. - :returns: A new UFL expression - """ - ifacet = self.interface.integral_type.startswith("interior_facet") - return preprocess_expression(expr, complex_mode=context.complex_mode, - do_apply_restrictions=ifacet) - - @property - def config(self): - config = {name: getattr(self.interface, name) - for name in ["ufl_cell", "index_cache", "scalar_type"]} - config["interface"] = self.interface - return config - - def cell_size(self): - return self.interface.cell_size(self.mt.restriction) - - def jacobian_at(self, point): - ps = PointSingleton(point) - expr = Jacobian(extract_unique_domain(self.mt.terminal)) - assert ps.expression.shape == (extract_unique_domain(expr).topological_dimension(), ) - if self.mt.restriction == '+': - expr = PositiveRestricted(expr) - elif self.mt.restriction == '-': - expr = NegativeRestricted(expr) - config = {"point_set": PointSingleton(point)} - config.update(self.config) - context = PointSetContext(**config) - expr = self.preprocess(expr, context) - return map_expr_dag(context.translator, expr) - - def detJ_at(self, point): - expr = JacobianDeterminant(extract_unique_domain(self.mt.terminal)) - if self.mt.restriction == '+': - expr = PositiveRestricted(expr) - elif self.mt.restriction == '-': - expr = NegativeRestricted(expr) - config = {"point_set": PointSingleton(point)} - config.update(self.config) - context = PointSetContext(**config) - expr = self.preprocess(expr, context) - return map_expr_dag(context.translator, expr) - - def reference_normals(self): - cell = self.interface.fiat_cell - sd = cell.get_spatial_dimension() - num_faces = len(cell.get_topology()[sd-1]) - - return gem.Literal(numpy.asarray([cell.compute_normal(i) for i in range(num_faces)])) - - def reference_edge_tangents(self): - cell = self.interface.fiat_cell - num_edges = len(cell.get_topology()[1]) - return gem.Literal(numpy.asarray([cell.compute_edge_tangent(i) for i in range(num_edges)])) - - def physical_tangents(self): - cell = self.interface.fiat_cell - sd = cell.get_spatial_dimension() - num_edges = len(cell.get_topology()[1]) - els = self.physical_edge_lengths() - rts = gem.ListTensor([cell.compute_tangents(1, i)[0] / els[i] for i in range(num_edges)]) - jac = self.jacobian_at(cell.make_points(sd, 0, sd+1)[0]) - - return rts @ jac.T - - def physical_normals(self): - cell = self.interface.fiat_cell - if not (isinstance(cell, UFCSimplex) and cell.get_dimension() == 2): - raise NotImplementedError("Can't do physical normals on that cell yet") - - num_edges = len(cell.get_topology()[1]) - pts = self.physical_tangents() - return gem.ListTensor([[pts[i, 1], -1*pts[i, 0]] for i in range(num_edges)]) - - def physical_edge_lengths(self): - expr = ufl.classes.CellEdgeVectors(extract_unique_domain(self.mt.terminal)) - if self.mt.restriction == '+': - expr = PositiveRestricted(expr) - elif self.mt.restriction == '-': - expr = NegativeRestricted(expr) - - cell = self.interface.fiat_cell - sd = cell.get_spatial_dimension() - num_edges = len(cell.get_topology()[1]) - expr = ufl.as_vector([ufl.sqrt(ufl.dot(expr[i, :], expr[i, :])) for i in range(num_edges)]) - config = {"point_set": PointSingleton(cell.make_points(sd, 0, sd+1)[0])} - config.update(self.config) - context = PointSetContext(**config) - expr = self.preprocess(expr, context) - 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(extract_unique_domain(self.mt.terminal)) - point_shape, = point_set.expression.shape - if entity is not None: - e, _ = entity - assert point_shape == e - else: - assert point_shape == extract_unique_domain(expr).topological_dimension() - if self.mt.restriction == '+': - expr = PositiveRestricted(expr) - elif self.mt.restriction == '-': - expr = NegativeRestricted(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) - expr = self.preprocess(expr, context) - 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), NeedsCoordinateMappingElement) - - -class PointSetContext(ContextBase): - """Context for compile-time known evaluation points.""" - - keywords = ContextBase.keywords + ( - 'quadrature_degree', - 'quadrature_rule', - 'point_set', - 'weight_expr', - ) - - @cached_property - def quadrature_rule(self): - integration_cell = self.fiat_cell.construct_subelement(self.integration_dim) - return make_quadrature(integration_cell, self.quadrature_degree) - - @cached_property - def point_set(self): - return self.quadrature_rule.point_set - - @cached_property - def point_indices(self): - return self.point_set.indices - - @cached_property - def point_expr(self): - return self.point_set.expression - - @cached_property - def weight_expr(self): - return self.quadrature_rule.weight_expression - - def basis_evaluation(self, finat_element, mt, entity_id): - return finat_element.basis_evaluation(mt.local_derivatives, - self.point_set, - (self.integration_dim, entity_id), - coordinate_mapping=CoordinateMapping(mt, self)) - - -class GemPointContext(ContextBase): - """Context for evaluation at arbitrary reference points.""" - - keywords = ContextBase.keywords + ( - 'point_indices', - 'point_expr', - 'weight_expr', - ) - - def basis_evaluation(self, finat_element, mt, entity_id): - return finat_element.point_evaluation(mt.local_derivatives, - self.point_expr, - (self.integration_dim, entity_id)) - - -class Translator(MultiFunction, ModifiedTerminalMixin, ufl2gem.Mixin): - """Multifunction for translating UFL -> GEM. Incorporates ufl2gem.Mixin, and - dispatches on terminal type when reaching modified terminals.""" - - def __init__(self, context): - # MultiFunction.__init__ does not call further __init__ - # methods, but ufl2gem.Mixin must be initialised. - # (ModifiedTerminalMixin requires no initialisation.) - MultiFunction.__init__(self) - ufl2gem.Mixin.__init__(self) - - # Need context during translation! - self.context = context - - # We just use the provided quadrature rule to - # perform the integration. - # Can't put these in the ufl2gem mixin, since they (unlike - # everything else) want access to the translation context. - def cell_avg(self, o): - if self.context.integral_type != "cell": - # Need to create a cell-based quadrature rule and - # translate the expression using that (c.f. CellVolume - # below). - raise NotImplementedError("CellAvg on non-cell integrals not yet implemented") - integrand, = o.ufl_operands - domain = extract_unique_domain(o) - measure = ufl.Measure(self.context.integral_type, domain=domain) - integrand, degree, argument_multiindices = entity_avg(integrand / CellVolume(domain), measure, self.context.argument_multiindices) - - config = {name: getattr(self.context, name) - for name in ["ufl_cell", "index_cache", "scalar_type"]} - config.update(quadrature_degree=degree, interface=self.context, - argument_multiindices=argument_multiindices) - expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True) - return expr - - def facet_avg(self, o): - if self.context.integral_type == "cell": - raise ValueError("Can't take FacetAvg in cell integral") - integrand, = o.ufl_operands - domain = extract_unique_domain(o) - measure = ufl.Measure(self.context.integral_type, domain=domain) - integrand, degree, argument_multiindices = entity_avg(integrand / FacetArea(domain), measure, self.context.argument_multiindices) - - config = {name: getattr(self.context, name) - for name in ["ufl_cell", "index_cache", "scalar_type", - "integration_dim", "entity_ids", - "integral_type"]} - config.update(quadrature_degree=degree, interface=self.context, - argument_multiindices=argument_multiindices) - expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True) - return expr - - def modified_terminal(self, o): - """Overrides the modified terminal handler from - :class:`ModifiedTerminalMixin`.""" - mt = analyse_modified_terminal(o) - return translate(mt.terminal, mt, self.context) - - -@singledispatch -def translate(terminal, mt, ctx): - """Translates modified terminals into GEM. - - :arg terminal: terminal, for dispatching - :arg mt: analysed modified terminal - :arg ctx: translator context - :returns: GEM translation of the modified terminal - """ - raise AssertionError("Cannot handle terminal type: %s" % type(terminal)) - - -@translate.register(QuadratureWeight) -def translate_quadratureweight(terminal, mt, ctx): - return ctx.weight_expr - - -@translate.register(GeometricQuantity) -def translate_geometricquantity(terminal, mt, ctx): - raise NotImplementedError("Cannot handle geometric quantity type: %s" % type(terminal)) - - -@translate.register(CellOrientation) -def translate_cell_orientation(terminal, mt, ctx): - return ctx.cell_orientation(mt.restriction) - - -@translate.register(ReferenceCellVolume) -def translate_reference_cell_volume(terminal, mt, ctx): - return gem.Literal(ctx.fiat_cell.volume()) - - -@translate.register(ReferenceFacetVolume) -def translate_reference_facet_volume(terminal, mt, ctx): - assert ctx.integral_type != "cell" - # Sum of quadrature weights is entity volume - return gem.optimise.aggressive_unroll(gem.index_sum(ctx.weight_expr, - ctx.point_indices)) - - -@translate.register(CellFacetJacobian) -def translate_cell_facet_jacobian(terminal, mt, ctx): - cell = ctx.fiat_cell - facet_dim = ctx.integration_dim - assert facet_dim != cell.get_dimension() - - def callback(entity_id): - return gem.Literal(make_cell_facet_jacobian(cell, facet_dim, entity_id)) - return ctx.entity_selector(callback, mt.restriction) - - -def make_cell_facet_jacobian(cell, facet_dim, facet_i): - facet_cell = cell.construct_subelement(facet_dim) - xs = facet_cell.get_vertices() - ys = cell.get_vertices_of_subcomplex(cell.get_topology()[facet_dim][facet_i]) - - # Use first 'dim' points to make an affine mapping - dim = cell.get_spatial_dimension() - A, b = make_affine_mapping(xs[:dim], ys[:dim]) - - for x, y in zip(xs[dim:], ys[dim:]): - # The rest of the points are checked to make sure the - # mapping really *is* affine. - assert numpy.allclose(y, A.dot(x) + b) - - return A - - -@translate.register(ReferenceNormal) -def translate_reference_normal(terminal, mt, ctx): - def callback(facet_i): - n = ctx.fiat_cell.compute_reference_normal(ctx.integration_dim, facet_i) - return gem.Literal(n) - return ctx.entity_selector(callback, mt.restriction) - - -@translate.register(ReferenceCellEdgeVectors) -def translate_reference_cell_edge_vectors(terminal, mt, ctx): - from FIAT.reference_element import \ - TensorProductCell as fiat_TensorProductCell - fiat_cell = ctx.fiat_cell - if isinstance(fiat_cell, fiat_TensorProductCell): - raise NotImplementedError("ReferenceCellEdgeVectors not implemented on TensorProductElements yet") - - nedges = len(fiat_cell.get_topology()[1]) - vecs = numpy.vstack(map(fiat_cell.compute_edge_tangent, range(nedges))) - assert vecs.shape == terminal.ufl_shape - return gem.Literal(vecs) - - -@translate.register(CellCoordinate) -def translate_cell_coordinate(terminal, mt, ctx): - if ctx.integration_dim == ctx.fiat_cell.get_dimension(): - return ctx.point_expr - - # This destroys the structure of the quadrature points, but since - # this code path is only used to implement CellCoordinate in facet - # integrals, hopefully it does not matter much. - ps = ctx.point_set - point_shape = tuple(index.extent for index in ps.indices) - - def callback(entity_id): - t = ctx.fiat_cell.get_entity_transform(ctx.integration_dim, entity_id) - data = numpy.asarray(list(map(t, ps.points))) - return gem.Literal(data.reshape(point_shape + data.shape[1:])) - - return gem.partial_indexed(ctx.entity_selector(callback, mt.restriction), - ps.indices) - - -@translate.register(FacetCoordinate) -def translate_facet_coordinate(terminal, mt, ctx): - assert ctx.integration_dim != ctx.fiat_cell.get_dimension() - return ctx.point_expr - - -@translate.register(SpatialCoordinate) -def translate_spatialcoordinate(terminal, mt, ctx): - # Replace terminal with a Coefficient - terminal = ctx.coordinate(extract_unique_domain(terminal)) - # Get back to reference space - terminal = preprocess_expression(terminal, complex_mode=ctx.complex_mode) - # Rebuild modified terminal - expr = construct_modified_terminal(mt, terminal) - # Translate replaced UFL snippet - return ctx.translator(expr) - - -class CellVolumeKernelInterface(ProxyKernelInterface): - # Since CellVolume is evaluated as a cell integral, we must ensure - # that the right restriction is applied when it is used in an - # interior facet integral. This proxy diverts coefficient - # translation to use a specified restriction. - - def __init__(self, wrapee, restriction): - ProxyKernelInterface.__init__(self, wrapee) - self.restriction = restriction - - def coefficient(self, ufl_coefficient, r): - assert r is None - return self._wrapee.coefficient(ufl_coefficient, self.restriction) - - -@translate.register(CellVolume) -def translate_cellvolume(terminal, mt, ctx): - integrand, degree = one_times(ufl.dx(domain=extract_unique_domain(terminal))) - interface = CellVolumeKernelInterface(ctx, mt.restriction) - - config = {name: getattr(ctx, name) - for name in ["ufl_cell", "index_cache", "scalar_type"]} - config.update(interface=interface, quadrature_degree=degree) - expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True) - return expr - - -@translate.register(FacetArea) -def translate_facetarea(terminal, mt, ctx): - assert ctx.integral_type != 'cell' - domain = extract_unique_domain(terminal) - integrand, degree = one_times(ufl.Measure(ctx.integral_type, domain=domain)) - - config = {name: getattr(ctx, name) - for name in ["ufl_cell", "integration_dim", "scalar_type", - "entity_ids", "index_cache"]} - config.update(interface=ctx, quadrature_degree=degree) - expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True) - return expr - - -@translate.register(CellOrigin) -def translate_cellorigin(terminal, mt, ctx): - domain = extract_unique_domain(terminal) - coords = SpatialCoordinate(domain) - expression = construct_modified_terminal(mt, coords) - point_set = PointSingleton((0.0,) * domain.topological_dimension()) - - config = {name: getattr(ctx, name) - for name in ["ufl_cell", "index_cache", "scalar_type"]} - config.update(interface=ctx, point_set=point_set) - context = PointSetContext(**config) - return context.translator(expression) - - -@translate.register(CellVertices) -def translate_cell_vertices(terminal, mt, ctx): - coords = SpatialCoordinate(extract_unique_domain(terminal)) - ufl_expr = construct_modified_terminal(mt, coords) - ps = PointSet(numpy.array(ctx.fiat_cell.get_vertices())) - - config = {name: getattr(ctx, name) - for name in ["ufl_cell", "index_cache", "scalar_type"]} - config.update(interface=ctx, point_set=ps) - context = PointSetContext(**config) - expr = context.translator(ufl_expr) - - # Wrap up point (vertex) index - c = gem.Index() - return gem.ComponentTensor(gem.Indexed(expr, (c,)), ps.indices + (c,)) - - -@translate.register(CellEdgeVectors) -def translate_cell_edge_vectors(terminal, mt, ctx): - # WARNING: Assumes straight edges! - coords = CellVertices(extract_unique_domain(terminal)) - ufl_expr = construct_modified_terminal(mt, coords) - cell_vertices = ctx.translator(ufl_expr) - - e = gem.Index() - c = gem.Index() - expr = gem.ListTensor([ - gem.Sum(gem.Indexed(cell_vertices, (u, c)), - gem.Product(gem.Literal(-1), - gem.Indexed(cell_vertices, (v, c)))) - for _, (u, v) in sorted(ctx.fiat_cell.get_topology()[1].items()) - ]) - return gem.ComponentTensor(gem.Indexed(expr, (e,)), (e, c)) - - -def fiat_to_ufl(fiat_dict, order): - # All derivative multiindices must be of the same dimension. - dimension, = set(len(alpha) for alpha in fiat_dict.keys()) - - # All derivative tables must have the same shape. - shape, = set(table.shape for table in fiat_dict.values()) - sigma = tuple(gem.Index(extent=extent) for extent in shape) - - # Convert from FIAT to UFL format - eye = numpy.eye(dimension, dtype=int) - tensor = numpy.empty((dimension,) * order, dtype=object) - for multiindex in numpy.ndindex(tensor.shape): - alpha = tuple(eye[multiindex, :].sum(axis=0)) - tensor[multiindex] = gem.Indexed(fiat_dict[alpha], sigma) - delta = tuple(gem.Index(extent=dimension) for _ in range(order)) - if order > 0: - tensor = gem.Indexed(gem.ListTensor(tensor), delta) - else: - tensor = tensor[()] - return gem.ComponentTensor(tensor, sigma + delta) - - -@translate.register(Argument) -def translate_argument(terminal, mt, ctx): - argument_multiindex = ctx.argument_multiindices[terminal.number()] - sigma = tuple(gem.Index(extent=d) for d in mt.expr.ufl_shape) - element = ctx.create_element(terminal.ufl_element(), restriction=mt.restriction) - - def callback(entity_id): - finat_dict = ctx.basis_evaluation(element, mt, entity_id) - # Filter out irrelevant derivatives - filtered_dict = {alpha: table - for alpha, table in finat_dict.items() - if sum(alpha) == mt.local_derivatives} - - # Change from FIAT to UFL arrangement - square = fiat_to_ufl(filtered_dict, mt.local_derivatives) - - # A numerical hack that FFC used to apply on FIAT tables still - # lives on after ditching FFC and switching to FInAT. - return ffc_rounding(square, ctx.epsilon) - table = ctx.entity_selector(callback, mt.restriction) - return gem.ComponentTensor(gem.Indexed(table, argument_multiindex + sigma), sigma) - - -@translate.register(TSFCConstantMixin) -def translate_constant_value(terminal, mt, ctx): - return ctx.constant(terminal) - - -@translate.register(Coefficient) -def translate_coefficient(terminal, mt, ctx): - vec = ctx.coefficient(terminal, mt.restriction) - - if terminal.ufl_element().family() == 'Real': - assert mt.local_derivatives == 0 - return vec - - element = ctx.create_element(terminal.ufl_element(), restriction=mt.restriction) - - # Collect FInAT tabulation for all entities - per_derivative = collections.defaultdict(list) - for entity_id in ctx.entity_ids: - finat_dict = ctx.basis_evaluation(element, mt, entity_id) - for alpha, table in finat_dict.items(): - # Filter out irrelevant derivatives - if sum(alpha) == mt.local_derivatives: - # A numerical hack that FFC used to apply on FIAT - # tables still lives on after ditching FFC and - # switching to FInAT. - table = ffc_rounding(table, ctx.epsilon) - per_derivative[alpha].append(table) - - # Merge entity tabulations for each derivative - if len(ctx.entity_ids) == 1: - def take_singleton(xs): - x, = xs # asserts singleton - return x - per_derivative = {alpha: take_singleton(tables) - for alpha, tables in per_derivative.items()} - else: - f = ctx.entity_number(mt.restriction) - per_derivative = {alpha: gem.select_expression(tables, f) - for alpha, tables in per_derivative.items()} - - # Coefficient evaluation - ctx.index_cache.setdefault(terminal.ufl_element(), element.get_indices()) - beta = ctx.index_cache[terminal.ufl_element()] - zeta = element.get_value_indices() - vec_beta, = gem.optimise.remove_componenttensors([gem.Indexed(vec, beta)]) - value_dict = {} - for alpha, table in per_derivative.items(): - table_qi = gem.Indexed(table, beta + zeta) - summands = [] - for var, expr in unconcatenate([(vec_beta, table_qi)], ctx.index_cache): - indices = tuple(i for i in var.index_ordering() if i not in ctx.unsummed_coefficient_indices) - value = gem.IndexSum(gem.Product(expr, var), indices) - summands.append(gem.optimise.contraction(value)) - optimised_value = gem.optimise.make_sum(summands) - value_dict[alpha] = gem.ComponentTensor(optimised_value, zeta) - - # Change from FIAT to UFL arrangement - result = fiat_to_ufl(value_dict, mt.local_derivatives) - assert result.shape == mt.expr.ufl_shape - assert set(result.free_indices) - ctx.unsummed_coefficient_indices <= set(ctx.point_indices) - - # Detect Jacobian of affine cells - if not result.free_indices and all(numpy.count_nonzero(node.array) <= 2 - for node in traversal((result,)) - if isinstance(node, gem.Literal)): - result = gem.optimise.aggressive_unroll(result) - return result - - -def compile_ufl(expression, context, interior_facet=False, point_sum=False): - """Translate a UFL expression to GEM. - - :arg expression: The UFL expression to compile. - :arg context: translation context - either a :class:`GemPointContext` - or :class:`PointSetContext` - :arg interior_facet: If ``true``, treat expression as an interior - facet integral (default ``False``) - :arg point_sum: If ``true``, return a `gem.IndexSum` of the final - gem expression along the ``context.point_indices`` (if present). - """ - - # Abs-simplification - expression = simplify_abs(expression, context.complex_mode) - if interior_facet: - expressions = [] - for rs in itertools.product(("+", "-"), repeat=len(context.argument_multiindices)): - expressions.append(map_expr_dag(PickRestriction(*rs), expression)) - else: - expressions = [expression] - - # Translate UFL to GEM, lowering finite element specific nodes - result = map_expr_dags(context.translator, expressions) - if point_sum: - result = [gem.index_sum(expr, context.point_indices) for expr in result] - return constant_fold_zero(result) diff --git a/tsfc/finatinterface.py b/tsfc/finatinterface.py deleted file mode 100644 index f9f461fc..00000000 --- a/tsfc/finatinterface.py +++ /dev/null @@ -1,359 +0,0 @@ -# This file was modified from FFC -# (http://bitbucket.org/fenics-project/ffc), copyright notice -# reproduced below. -# -# Copyright (C) 2009-2013 Kristian B. Oelgaard and Anders Logg -# -# This file is part of FFC. -# -# FFC 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. -# -# FFC 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 FFC. If not, see . - -import weakref -from functools import partial, singledispatch - -import FIAT -import finat -import finat.ufl -import ufl - -__all__ = ("as_fiat_cell", "create_base_element", - "create_element", "supported_elements") - - -supported_elements = { - # These all map directly to FInAT elements - "Bernstein": finat.Bernstein, - "Bernardi-Raugel": finat.BernardiRaugel, - "Bernardi-Raugel Bubble": finat.BernardiRaugelBubble, - "Brezzi-Douglas-Marini": finat.BrezziDouglasMarini, - "Brezzi-Douglas-Fortin-Marini": finat.BrezziDouglasFortinMarini, - "Bubble": finat.Bubble, - "FacetBubble": finat.FacetBubble, - "Crouzeix-Raviart": finat.CrouzeixRaviart, - "Discontinuous Lagrange": finat.DiscontinuousLagrange, - "Discontinuous Raviart-Thomas": lambda c, d: finat.DiscontinuousElement(finat.RaviartThomas(c, d)), - "Discontinuous Taylor": finat.DiscontinuousTaylor, - "Gauss-Legendre": finat.GaussLegendre, - "Gauss-Lobatto-Legendre": finat.GaussLobattoLegendre, - "HDiv Trace": finat.HDivTrace, - "Hellan-Herrmann-Johnson": finat.HellanHerrmannJohnson, - "Johnson-Mercier": finat.JohnsonMercier, - "Nonconforming Arnold-Winther": finat.ArnoldWintherNC, - "Conforming Arnold-Winther": finat.ArnoldWinther, - "Hermite": finat.Hermite, - "Kong-Mulder-Veldhuizen": finat.KongMulderVeldhuizen, - "Argyris": finat.Argyris, - "Hsieh-Clough-Tocher": finat.HsiehCloughTocher, - "QuadraticPowellSabin6": finat.QuadraticPowellSabin6, - "QuadraticPowellSabin12": finat.QuadraticPowellSabin12, - "Reduced-Hsieh-Clough-Tocher": finat.ReducedHsiehCloughTocher, - "Mardal-Tai-Winther": finat.MardalTaiWinther, - "Alfeld-Sorokina": finat.AlfeldSorokina, - "Arnold-Qin": finat.ArnoldQin, - "Reduced-Arnold-Qin": finat.ReducedArnoldQin, - "Christiansen-Hu": finat.ChristiansenHu, - "Guzman-Neilan 1st kind H1": finat.GuzmanNeilanFirstKindH1, - "Guzman-Neilan 2nd kind H1": finat.GuzmanNeilanSecondKindH1, - "Guzman-Neilan Bubble": finat.GuzmanNeilanBubble, - "Guzman-Neilan H1(div)": finat.GuzmanNeilanH1div, - "Morley": finat.Morley, - "Bell": finat.Bell, - "Lagrange": finat.Lagrange, - "Nedelec 1st kind H(curl)": finat.Nedelec, - "Nedelec 2nd kind H(curl)": finat.NedelecSecondKind, - "Raviart-Thomas": finat.RaviartThomas, - "Regge": finat.Regge, - "BDMCE": finat.BrezziDouglasMariniCubeEdge, - "BDMCF": finat.BrezziDouglasMariniCubeFace, - # These require special treatment below - "DQ": None, - "Q": None, - "RTCE": None, - "RTCF": None, - "NCE": None, - "NCF": None, - "Real": finat.Real, - "DPC": finat.DPC, - "S": finat.Serendipity, - "SminusF": finat.TrimmedSerendipityFace, - "SminusDiv": finat.TrimmedSerendipityDiv, - "SminusE": finat.TrimmedSerendipityEdge, - "SminusCurl": finat.TrimmedSerendipityCurl, - "DPC L2": finat.DPC, - "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 -element is supported, but must be handled specially because it doesn't -have a direct FInAT equivalent.""" - - -def as_fiat_cell(cell): - """Convert a ufl cell to a FIAT cell. - - :arg cell: the :class:`ufl.Cell` to convert.""" - if not isinstance(cell, ufl.AbstractCell): - raise ValueError("Expecting a UFL Cell") - return FIAT.ufc_cell(cell) - - -@singledispatch -def convert(element, **kwargs): - """Handler for converting UFL elements to FInAT elements. - - :arg element: The UFL element to convert. - - Do not use this function directly, instead call - :func:`create_element`.""" - if element.family() in supported_elements: - raise ValueError("Element %s supported, but no handler provided" % element) - raise ValueError("Unsupported element type %s" % type(element)) - - -# Base finite elements first -@convert.register(finat.ufl.FiniteElement) -def convert_finiteelement(element, **kwargs): - cell = as_fiat_cell(element.cell) - if element.family() == "Quadrature": - degree = element.degree() - scheme = element.quadrature_scheme() - if degree is None or scheme is None: - raise ValueError("Quadrature scheme and degree must be specified!") - - return finat.make_quadrature_element(cell, degree, scheme), set() - lmbda = supported_elements[element.family()] - if element.family() == "Real" and element.cell.cellname() in {"quadrilateral", "hexahedron"}: - lmbda = None - element = finat.ufl.FiniteElement("DQ", element.cell, 0) - if lmbda is None: - if element.cell.cellname() == "quadrilateral": - # Handle quadrilateral short names like RTCF and RTCE. - element = element.reconstruct(cell=quadrilateral_tpc) - elif element.cell.cellname() == "hexahedron": - # Handle hexahedron short names like NCF and NCE. - element = element.reconstruct(cell=hexahedron_tpc) - else: - raise ValueError("%s is supported, but handled incorrectly" % - element.family()) - finat_elem, deps = _create_element(element, **kwargs) - return finat.FlattenedDimensions(finat_elem), deps - - kind = element.variant() - if kind is None: - kind = 'spectral' # default variant - is_interval = element.cell.cellname() == 'interval' - - if element.family() in {"Raviart-Thomas", "Nedelec 1st kind H(curl)", - "Brezzi-Douglas-Marini", "Nedelec 2nd kind H(curl)", - "Argyris"}: - lmbda = partial(lmbda, variant=element.variant()) - elif element.family() == "Lagrange": - if kind == 'spectral': - lmbda = finat.GaussLobattoLegendre - elif kind.startswith('integral'): - lmbda = partial(finat.IntegratedLegendre, variant=kind) - elif kind in ['fdm', 'fdm_ipdg'] and is_interval: - lmbda = finat.FDMLagrange - elif kind == 'fdm_quadrature' and is_interval: - lmbda = finat.FDMQuadrature - elif kind == 'fdm_broken' and is_interval: - lmbda = finat.FDMBrokenH1 - elif kind == 'fdm_hermite' and is_interval: - lmbda = finat.FDMHermite - elif kind in ['demkowicz', 'fdm']: - lmbda = partial(finat.IntegratedLegendre, variant=kind) - elif kind in ['mgd', 'feec', 'qb', 'mse']: - degree = element.degree() - shift_axes = kwargs["shift_axes"] - restriction = kwargs["restriction"] - deps = {"shift_axes", "restriction"} - return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction), deps - else: - # Let FIAT handle the general case - lmbda = partial(finat.Lagrange, variant=kind) - elif element.family() in ["Discontinuous Lagrange", "Discontinuous Lagrange L2"]: - if kind == 'spectral': - lmbda = finat.GaussLegendre - elif kind.startswith('integral'): - lmbda = partial(finat.Legendre, variant=kind) - elif kind in ['fdm', 'fdm_quadrature'] and is_interval: - lmbda = finat.FDMDiscontinuousLagrange - elif kind == 'fdm_ipdg' and is_interval: - lmbda = lambda *args: finat.DiscontinuousElement(finat.FDMLagrange(*args)) - elif kind in 'fdm_broken' and is_interval: - lmbda = finat.FDMBrokenL2 - elif kind in ['demkowicz', 'fdm']: - lmbda = partial(finat.Legendre, variant=kind) - elif kind in ['mgd', 'feec', 'qb', 'mse']: - degree = element.degree() - shift_axes = kwargs["shift_axes"] - restriction = kwargs["restriction"] - deps = {"shift_axes", "restriction"} - return finat.RuntimeTabulated(cell, degree, variant=kind, shift_axes=shift_axes, restriction=restriction, continuous=False), deps - else: - # Let FIAT handle the general case - lmbda = partial(finat.DiscontinuousLagrange, variant=kind) - elif element.family() == ["DPC", "DPC L2", "S"]: - dim = element.cell.geometric_dimension() - if dim > 1: - element = element.reconstruct(cell=ufl.cell.hypercube(dim)) - - return lmbda(cell, element.degree()), set() - - -# Element modifiers and compound element types -@convert.register(finat.ufl.BrokenElement) -def convert_brokenelement(element, **kwargs): - finat_elem, deps = _create_element(element._element, **kwargs) - return finat.DiscontinuousElement(finat_elem), deps - - -@convert.register(finat.ufl.EnrichedElement) -def convert_enrichedelement(element, **kwargs): - elements, deps = zip(*[_create_element(elem, **kwargs) - for elem in element._elements]) - return finat.EnrichedElement(elements), set.union(*deps) - - -@convert.register(finat.ufl.NodalEnrichedElement) -def convert_nodalenrichedelement(element, **kwargs): - elements, deps = zip(*[_create_element(elem, **kwargs) - for elem in element._elements]) - return finat.NodalEnrichedElement(elements), set.union(*deps) - - -@convert.register(finat.ufl.MixedElement) -def convert_mixedelement(element, **kwargs): - elements, deps = zip(*[_create_element(elem, **kwargs) - for elem in element.sub_elements]) - return finat.MixedElement(elements), set.union(*deps) - - -@convert.register(finat.ufl.VectorElement) -@convert.register(finat.ufl.TensorElement) -def convert_tensorelement(element, **kwargs): - inner_elem, deps = _create_element(element.sub_elements[0], **kwargs) - shape = element.reference_value_shape - shape = shape[:len(shape) - len(inner_elem.value_shape)] - shape_innermost = kwargs["shape_innermost"] - return (finat.TensorFiniteElement(inner_elem, shape, not shape_innermost), - deps | {"shape_innermost"}) - - -@convert.register(finat.ufl.TensorProductElement) -def convert_tensorproductelement(element, **kwargs): - cell = element.cell - if type(cell) is not ufl.TensorProductCell: - raise ValueError("TensorProductElement not on TensorProductCell?") - shift_axes = kwargs["shift_axes"] - dim_offset = 0 - elements = [] - deps = set() - for elem in element.sub_elements: - kwargs["shift_axes"] = shift_axes + dim_offset - dim_offset += elem.cell.topological_dimension() - finat_elem, ds = _create_element(elem, **kwargs) - elements.append(finat_elem) - deps.update(ds) - return finat.TensorProductElement(elements), deps - - -@convert.register(finat.ufl.HDivElement) -def convert_hdivelement(element, **kwargs): - finat_elem, deps = _create_element(element._element, **kwargs) - return finat.HDivElement(finat_elem), deps - - -@convert.register(finat.ufl.HCurlElement) -def convert_hcurlelement(element, **kwargs): - finat_elem, deps = _create_element(element._element, **kwargs) - return finat.HCurlElement(finat_elem), deps - - -@convert.register(finat.ufl.WithMapping) -def convert_withmapping(element, **kwargs): - return _create_element(element.wrapee, **kwargs) - - -@convert.register(finat.ufl.RestrictedElement) -def convert_restrictedelement(element, **kwargs): - finat_elem, deps = _create_element(element._element, **kwargs) - return finat.RestrictedElement(finat_elem, element.restriction_domain()), deps - - -hexahedron_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval, ufl.interval) -quadrilateral_tpc = ufl.TensorProductCell(ufl.interval, ufl.interval) -_cache = weakref.WeakKeyDictionary() - - -def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction=None): - """Create a FInAT element (suitable for tabulating with) given a UFL element. - - :arg ufl_element: The UFL element to create a FInAT element from. - :arg shape_innermost: Vector/tensor indices come after basis function indices - :arg restriction: cell restriction in interior facet integrals - (only for runtime tabulated elements) - """ - finat_element, deps = _create_element(ufl_element, - shape_innermost=shape_innermost, - shift_axes=shift_axes, - restriction=restriction) - return finat_element - - -def _create_element(ufl_element, **kwargs): - """A caching wrapper around :py:func:`convert`. - - Takes a UFL element and an unspecified set of parameter options, - and returns the converted element with the set of keyword names - that were relevant for conversion. - """ - # Look up conversion in cache - try: - cache = _cache[ufl_element] - except KeyError: - _cache[ufl_element] = {} - cache = _cache[ufl_element] - - for key, finat_element in cache.items(): - # Cache hit if all relevant parameter values match. - if all(kwargs[param] == value for param, value in key): - return finat_element, set(param for param, value in key) - - # Convert if cache miss - if ufl_element.cell is None: - raise ValueError("Don't know how to build element when cell is not given") - - finat_element, deps = convert(ufl_element, **kwargs) - - # Store conversion in cache - key = frozenset((param, kwargs[param]) for param in deps) - cache[key] = finat_element - - # Forward result - return finat_element, deps - - -def create_base_element(ufl_element, **kwargs): - """Create a "scalar" base FInAT element given a UFL element. - Takes a UFL element and an unspecified set of parameter options, - and returns the converted element. - """ - finat_element = create_element(ufl_element, **kwargs) - if isinstance(finat_element, finat.TensorFiniteElement): - finat_element = finat_element.base_element - return finat_element diff --git a/tsfc/kernel_args.py b/tsfc/kernel_args.py deleted file mode 100644 index 1c79bdae..00000000 --- a/tsfc/kernel_args.py +++ /dev/null @@ -1,54 +0,0 @@ -import abc - - -class KernelArg(abc.ABC): - """Abstract base class wrapping a loopy argument. - - Defining this type system allows Firedrake (or other libraries) to - prepare arguments for the kernel without needing to worry about their - ordering. Instead this can be offloaded to tools such as - :func:`functools.singledispatch`. - """ - - def __init__(self, arg): - self.loopy_arg = arg - - @property - def dtype(self): - return self.loopy_arg.dtype - - -class OutputKernelArg(KernelArg): - ... - - -class CoordinatesKernelArg(KernelArg): - ... - - -class CoefficientKernelArg(KernelArg): - ... - - -class ConstantKernelArg(KernelArg): - ... - - -class CellOrientationsKernelArg(KernelArg): - ... - - -class CellSizesKernelArg(KernelArg): - ... - - -class TabulationKernelArg(KernelArg): - ... - - -class ExteriorFacetKernelArg(KernelArg): - ... - - -class InteriorFacetKernelArg(KernelArg): - ... diff --git a/tsfc/kernel_interface/__init__.py b/tsfc/kernel_interface/__init__.py deleted file mode 100644 index fc194260..00000000 --- a/tsfc/kernel_interface/__init__.py +++ /dev/null @@ -1,47 +0,0 @@ -from abc import ABCMeta, abstractmethod, abstractproperty - -from gem.utils import make_proxy_class - - -class KernelInterface(metaclass=ABCMeta): - """Abstract interface for accessing the GEM expressions corresponding - to kernel arguments.""" - - @abstractmethod - def coordinate(self, ufl_domain): - """A function that maps :class:`ufl.Domain`s to coordinate - :class:`ufl.Coefficient`s.""" - - @abstractmethod - def coefficient(self, ufl_coefficient, restriction): - """A function that maps :class:`ufl.Coefficient`s to GEM - expressions.""" - - @abstractmethod - def constant(self, const): - """Return the GEM expression corresponding to the constant.""" - - @abstractmethod - def cell_orientation(self, restriction): - """Cell orientation as a GEM expression.""" - - @abstractmethod - def cell_size(self, restriction): - """Mesh cell size as a GEM expression. Shape (nvertex, ) in FIAT vertex ordering.""" - - @abstractmethod - def entity_number(self, restriction): - """Facet or vertex number as a GEM index.""" - - @abstractmethod - def create_element(self, element, **kwargs): - """Create a FInAT element (suitable for tabulating with) given - a UFL element.""" - - @abstractproperty - def unsummed_coefficient_indices(self): - """A set of indices that coefficient evaluation should not sum over. - Used for macro-cell integration.""" - - -ProxyKernelInterface = make_proxy_class('ProxyKernelInterface', KernelInterface) diff --git a/tsfc/kernel_interface/common.py b/tsfc/kernel_interface/common.py deleted file mode 100644 index 18fd363f..00000000 --- a/tsfc/kernel_interface/common.py +++ /dev/null @@ -1,540 +0,0 @@ -import collections -import operator -import string -from functools import reduce -from itertools import chain, product - -import gem -import gem.impero_utils as impero_utils -import numpy -from FIAT.reference_element import TensorProductCell -from finat.cell_tools import max_complex -from finat.quadrature import AbstractQuadratureRule, make_quadrature -from gem.node import traversal -from gem.optimise import constant_fold_zero -from gem.optimise import remove_componenttensors as prune -from gem.utils import cached_property -from numpy import asarray -from tsfc import fem, ufl_utils -from tsfc.finatinterface import as_fiat_cell, create_element -from tsfc.kernel_interface import KernelInterface -from tsfc.logging import logger -from ufl.utils.sequences import max_degree - - -class KernelBuilderBase(KernelInterface): - """Helper class for building local assembly kernels.""" - - def __init__(self, scalar_type, interior_facet=False): - """Initialise a kernel builder. - - :arg interior_facet: kernel accesses two cells - """ - assert isinstance(interior_facet, bool) - self.scalar_type = scalar_type - self.interior_facet = interior_facet - - self.prepare = [] - self.finalise = [] - - # Coordinates - self.domain_coordinate = {} - - # Coefficients - self.coefficient_map = collections.OrderedDict() - - # Constants - self.constant_map = collections.OrderedDict() - - @cached_property - def unsummed_coefficient_indices(self): - return frozenset() - - def coordinate(self, domain): - return self.domain_coordinate[domain] - - def coefficient(self, ufl_coefficient, restriction): - """A function that maps :class:`ufl.Coefficient`s to GEM - expressions.""" - kernel_arg = self.coefficient_map[ufl_coefficient] - if ufl_coefficient.ufl_element().family() == 'Real': - return kernel_arg - elif not self.interior_facet: - return kernel_arg - else: - return kernel_arg[{'+': 0, '-': 1}[restriction]] - - def constant(self, const): - return self.constant_map[const] - - def cell_orientation(self, restriction): - """Cell orientation as a GEM expression.""" - f = {None: 0, '+': 0, '-': 1}[restriction] - # Assume self._cell_orientations tuple is set up at this point. - co_int = self._cell_orientations[f] - return gem.Conditional(gem.Comparison("==", co_int, gem.Literal(1)), - gem.Literal(-1), - gem.Conditional(gem.Comparison("==", co_int, gem.Zero()), - gem.Literal(1), - gem.Literal(numpy.nan))) - - def cell_size(self, restriction): - if not hasattr(self, "_cell_sizes"): - raise RuntimeError("Haven't called set_cell_sizes") - if self.interior_facet: - return self._cell_sizes[{'+': 0, '-': 1}[restriction]] - else: - return self._cell_sizes - - def entity_number(self, restriction): - """Facet or vertex number as a GEM index.""" - # Assume self._entity_number dict is set up at this point. - return self._entity_number[restriction] - - def apply_glue(self, prepare=None, finalise=None): - """Append glue code for operations that are not handled in the - GEM abstraction. - - Current uses: mixed interior facet mess - - :arg prepare: code snippets to be prepended to the kernel - :arg finalise: code snippets to be appended to the kernel - """ - if prepare is not None: - self.prepare.extend(prepare) - if finalise is not None: - self.finalise.extend(finalise) - - def register_requirements(self, ir): - """Inspect what is referenced by the IR that needs to be - provided by the kernel interface. - - :arg ir: multi-root GEM expression DAG - """ - # Nothing is required by default - pass - - -class KernelBuilderMixin(object): - """Mixin for KernelBuilder classes.""" - - def compile_integrand(self, integrand, params, ctx): - """Compile UFL integrand. - - :arg integrand: UFL integrand. - :arg params: a dict containing "quadrature_rule". - :arg ctx: context created with :meth:`create_context` method. - - See :meth:`create_context` for typical calling sequence. - """ - # Split Coefficients - if self.coefficient_split: - integrand = ufl_utils.split_coefficients(integrand, self.coefficient_split) - # Compile: ufl -> gem - info = self.integral_data_info - functions = list(info.arguments) + [self.coordinate(info.domain)] + list(info.coefficients) - set_quad_rule(params, info.domain.ufl_cell(), info.integral_type, functions) - quad_rule = params["quadrature_rule"] - config = self.fem_config() - config['argument_multiindices'] = self.argument_multiindices - config['quadrature_rule'] = quad_rule - config['index_cache'] = ctx['index_cache'] - expressions = fem.compile_ufl(integrand, - fem.PointSetContext(**config), - interior_facet=self.interior_facet) - ctx['quadrature_indices'].extend(quad_rule.point_set.indices) - return expressions - - def construct_integrals(self, integrand_expressions, params): - """Construct integrals from integrand expressions. - - :arg integrand_expressions: gem expressions for integrands. - :arg params: a dict containing "mode" and "quadrature_rule". - - integrand_expressions must be indexed with :attr:`argument_multiindices`; - these gem expressions are obtained by calling :meth:`compile_integrand` - method or by modifying the gem expressions returned by - :meth:`compile_integrand`. - - See :meth:`create_context` for typical calling sequence. - """ - mode = pick_mode(params["mode"]) - return mode.Integrals(integrand_expressions, - params["quadrature_rule"].point_set.indices, - self.argument_multiindices, - params) - - def stash_integrals(self, reps, params, ctx): - """Stash integral representations in ctx. - - :arg reps: integral representations. - :arg params: a dict containing "mode". - :arg ctx: context in which reps are stored. - - See :meth:`create_context` for typical calling sequence. - """ - mode = pick_mode(params["mode"]) - mode_irs = ctx['mode_irs'] - mode_irs.setdefault(mode, collections.OrderedDict()) - for var, rep in zip(self.return_variables, reps): - mode_irs[mode].setdefault(var, []).append(rep) - - def compile_gem(self, ctx): - """Compile gem representation of integrals to impero_c. - - :arg ctx: the context containing the gem representation of integrals. - :returns: a tuple of impero_c, oriented, needs_cell_sizes, tabulations - required to finally construct a kernel in :meth:`construct_kernel`. - - See :meth:`create_context` for typical calling sequence. - """ - # Finalise mode representations into a set of assignments - mode_irs = ctx['mode_irs'] - - assignments = [] - for mode, var_reps in mode_irs.items(): - assignments.extend(mode.flatten(var_reps.items(), ctx['index_cache'])) - - if assignments: - return_variables, expressions = zip(*assignments) - else: - return_variables = [] - expressions = [] - expressions = constant_fold_zero(expressions) - - # Need optimised roots - options = dict(reduce(operator.and_, - [mode.finalise_options.items() - for mode in mode_irs.keys()])) - expressions = impero_utils.preprocess_gem(expressions, **options) - - # Let the kernel interface inspect the optimised IR to register - # what kind of external data is required (e.g., cell orientations, - # cell sizes, etc.). - oriented, needs_cell_sizes, tabulations = self.register_requirements(expressions) - - # Extract Variables that are actually used - active_variables = gem.extract_type(expressions, gem.Variable) - # Construct ImperoC - assignments = list(zip(return_variables, expressions)) - index_ordering = get_index_ordering(ctx['quadrature_indices'], return_variables) - try: - impero_c = impero_utils.compile_gem(assignments, index_ordering, remove_zeros=True) - except impero_utils.NoopError: - impero_c = None - return impero_c, oriented, needs_cell_sizes, tabulations, active_variables - - def fem_config(self): - """Return a dictionary used with fem.compile_ufl. - - One needs to update this dictionary with "argument_multiindices", - "quadrature_rule", and "index_cache" before using this with - fem.compile_ufl. - """ - info = self.integral_data_info - integral_type = info.integral_type - cell = info.domain.ufl_cell() - fiat_cell = as_fiat_cell(cell) - integration_dim, entity_ids = lower_integral_type(fiat_cell, integral_type) - return dict(interface=self, - ufl_cell=cell, - integral_type=integral_type, - integration_dim=integration_dim, - entity_ids=entity_ids, - scalar_type=self.fem_scalar_type) - - def create_context(self): - """Create builder context. - - *index_cache* - - Map from UFL FiniteElement objects to multiindices. - This is so we reuse Index instances when evaluating the same - coefficient multiple times with the same table. - - We also use the same dict for the unconcatenate index cache, - which maps index objects to tuples of multiindices. These two - caches shall never conflict as their keys have different types - (UFL finite elements vs. GEM index objects). - - *quadrature_indices* - - List of quadrature indices used. - - *mode_irs* - - Dict for mode representations. - - For each set of integrals to make a kernel for (i,e., - `integral_data.integrals`), one must first create a ctx object by - calling :meth:`create_context` method. - This ctx object collects objects associated with the integrals that - are eventually used to construct the kernel. - The following is a typical calling sequence: - - .. code-block:: python3 - - builder = KernelBuilder(...) - params = {"mode": "spectral"} - ctx = builder.create_context() - for integral in integral_data.integrals: - integrand = integral.integrand() - integrand_exprs = builder.compile_integrand(integrand, params, ctx) - integral_exprs = builder.construct_integrals(integrand_exprs, params) - builder.stash_integrals(integral_exprs, params, ctx) - kernel = builder.construct_kernel(kernel_name, ctx) - - """ - return {'index_cache': {}, - 'quadrature_indices': [], - 'mode_irs': collections.OrderedDict()} - - -def set_quad_rule(params, cell, integral_type, functions): - # Check if the integral has a quad degree attached, otherwise use - # the estimated polynomial degree attached by compute_form_data - try: - quadrature_degree = params["quadrature_degree"] - except KeyError: - quadrature_degree = params["estimated_polynomial_degree"] - function_degrees = [f.ufl_function_space().ufl_element().degree() - for f in functions] - if all((asarray(quadrature_degree) > 10 * asarray(degree)).all() - for degree in function_degrees): - logger.warning("Estimated quadrature degree %s more " - "than tenfold greater than any " - "argument/coefficient degree (max %s)", - quadrature_degree, max_degree(function_degrees)) - quad_rule = params.get("quadrature_rule", "default") - if isinstance(quad_rule, str): - scheme = quad_rule - fiat_cell = as_fiat_cell(cell) - finat_elements = set(create_element(f.ufl_element()) for f in functions - if f.ufl_element().family() != "Real") - fiat_cells = [fiat_cell] + [finat_el.complex for finat_el in finat_elements] - fiat_cell = max_complex(fiat_cells) - - integration_dim, _ = lower_integral_type(fiat_cell, integral_type) - integration_cell = fiat_cell.construct_subcomplex(integration_dim) - quad_rule = make_quadrature(integration_cell, quadrature_degree, scheme=scheme) - params["quadrature_rule"] = quad_rule - - if not isinstance(quad_rule, AbstractQuadratureRule): - raise ValueError("Expected to find a QuadratureRule object, not a %s" % - type(quad_rule)) - - -def get_index_ordering(quadrature_indices, return_variables): - split_argument_indices = tuple(chain(*[var.index_ordering() - for var in return_variables])) - return tuple(quadrature_indices) + split_argument_indices - - -def get_index_names(quadrature_indices, argument_multiindices, index_cache): - index_names = [] - - def name_index(index, name): - index_names.append((index, name)) - if index in index_cache: - for multiindex, suffix in zip(index_cache[index], - string.ascii_lowercase): - name_multiindex(multiindex, name + suffix) - - def name_multiindex(multiindex, name): - if len(multiindex) == 1: - name_index(multiindex[0], name) - else: - for i, index in enumerate(multiindex): - name_index(index, name + str(i)) - - name_multiindex(quadrature_indices, 'ip') - for multiindex, name in zip(argument_multiindices, ['j', 'k']): - name_multiindex(multiindex, name) - return index_names - - -def lower_integral_type(fiat_cell, integral_type): - """Lower integral type into the dimension of the integration - subentity and a list of entity numbers for that dimension. - - :arg fiat_cell: FIAT reference cell - :arg integral_type: integral type (string) - """ - vert_facet_types = ['exterior_facet_vert', 'interior_facet_vert'] - horiz_facet_types = ['exterior_facet_bottom', 'exterior_facet_top', 'interior_facet_horiz'] - - dim = fiat_cell.get_dimension() - if integral_type == 'cell': - integration_dim = dim - elif integral_type in ['exterior_facet', 'interior_facet']: - if isinstance(fiat_cell, TensorProductCell): - raise ValueError("{} integral cannot be used with a TensorProductCell; need to distinguish between vertical and horizontal contributions.".format(integral_type)) - integration_dim = dim - 1 - elif integral_type == 'vertex': - integration_dim = 0 - elif integral_type in vert_facet_types + horiz_facet_types: - # Extrusion case - if not isinstance(fiat_cell, TensorProductCell): - raise ValueError("{} integral requires a TensorProductCell.".format(integral_type)) - basedim, extrdim = dim - assert extrdim == 1 - - if integral_type in vert_facet_types: - integration_dim = (basedim - 1, 1) - elif integral_type in horiz_facet_types: - integration_dim = (basedim, 0) - else: - raise NotImplementedError("integral type %s not supported" % integral_type) - - if integral_type == 'exterior_facet_bottom': - entity_ids = [0] - elif integral_type == 'exterior_facet_top': - entity_ids = [1] - else: - entity_ids = list(range(len(fiat_cell.get_topology()[integration_dim]))) - - return integration_dim, entity_ids - - -def pick_mode(mode): - "Return one of the specialized optimisation modules from a mode string." - try: - from firedrake_citations import Citations - cites = {"vanilla": ("Homolya2017", ), - "coffee": ("Luporini2016", "Homolya2017", ), - "spectral": ("Luporini2016", "Homolya2017", "Homolya2017a"), - "tensor": ("Kirby2006", "Homolya2017", )} - for c in cites[mode]: - Citations().register(c) - except ImportError: - pass - if mode == "vanilla": - import tsfc.vanilla as m - elif mode == "coffee": - import tsfc.coffee_mode as m - elif mode == "spectral": - import tsfc.spectral as m - elif mode == "tensor": - import tsfc.tensor as m - else: - raise ValueError("Unknown mode: {}".format(mode)) - return m - - -def check_requirements(ir): - """Look for cell orientations, cell sizes, and collect tabulations - in one pass.""" - cell_orientations = False - cell_sizes = False - rt_tabs = {} - for node in traversal(ir): - if isinstance(node, gem.Variable): - if node.name == "cell_orientations": - cell_orientations = True - elif node.name == "cell_sizes": - cell_sizes = True - elif node.name.startswith("rt_"): - rt_tabs[node.name] = node.shape - return cell_orientations, cell_sizes, tuple(sorted(rt_tabs.items())) - - -def prepare_constant(constant, number): - """Bridges the kernel interface and the GEM abstraction for - Constants. - - :arg constant: Firedrake Constant - :arg number: Value to uniquely identify the constant - :returns: (funarg, expression) - expression - GEM expression referring to the Constant value(s) - """ - value_size = numpy.prod(constant.ufl_shape, dtype=int) - return gem.reshape(gem.Variable(f"c_{number}", (value_size,)), - constant.ufl_shape) - - -def prepare_coefficient(coefficient, name, interior_facet=False): - """Bridges the kernel interface and the GEM abstraction for - Coefficients. - - :arg coefficient: UFL Coefficient - :arg name: unique name to refer to the Coefficient in the kernel - :arg interior_facet: interior facet integral? - :returns: (funarg, expression) - expression - GEM expression referring to the Coefficient - values - """ - assert isinstance(interior_facet, bool) - - if coefficient.ufl_element().family() == 'Real': - # Constant - value_size = coefficient.ufl_element().value_size - expression = gem.reshape(gem.Variable(name, (value_size,)), - coefficient.ufl_shape) - return expression - - finat_element = create_element(coefficient.ufl_element()) - shape = finat_element.index_shape - size = numpy.prod(shape, dtype=int) - - if not interior_facet: - expression = gem.reshape(gem.Variable(name, (size,)), shape) - else: - varexp = gem.Variable(name, (2 * size,)) - plus = gem.view(varexp, slice(size)) - minus = gem.view(varexp, slice(size, 2 * size)) - expression = (gem.reshape(plus, shape), gem.reshape(minus, shape)) - return expression - - -def prepare_arguments(arguments, multiindices, interior_facet=False, diagonal=False): - """Bridges the kernel interface and the GEM abstraction for - Arguments. Vector Arguments are rearranged here for interior - facet integrals. - - :arg arguments: UFL Arguments - :arg multiindices: Argument multiindices - :arg interior_facet: interior facet integral? - :arg diagonal: Are we assembling the diagonal of a rank-2 element tensor? - :returns: (funarg, expression) - expressions - GEM expressions referring to the argument - tensor - """ - assert isinstance(interior_facet, bool) - - if len(arguments) == 0: - # No arguments - expression = gem.Indexed(gem.Variable("A", (1,)), (0,)) - return (expression, ) - - elements = tuple(create_element(arg.ufl_element()) for arg in arguments) - shapes = tuple(element.index_shape for element in elements) - - if diagonal: - if len(arguments) != 2: - raise ValueError("Diagonal only for 2-forms") - try: - element, = set(elements) - except ValueError: - raise ValueError("Diagonal only for diagonal blocks (test and trial spaces the same)") - - elements = (element, ) - shapes = tuple(element.index_shape for element in elements) - multiindices = multiindices[:1] - - def expression(restricted): - return gem.Indexed(gem.reshape(restricted, *shapes), - tuple(chain(*multiindices))) - - u_shape = numpy.array([numpy.prod(shape, dtype=int) for shape in shapes]) - if interior_facet: - c_shape = tuple(2 * u_shape) - slicez = [[slice(r * s, (r + 1) * s) - for r, s in zip(restrictions, u_shape)] - for restrictions in product((0, 1), repeat=len(arguments))] - else: - c_shape = tuple(u_shape) - slicez = [[slice(s) for s in u_shape]] - - varexp = gem.Variable("A", c_shape) - expressions = [expression(gem.view(varexp, *slices)) for slices in slicez] - return tuple(prune(expressions)) diff --git a/tsfc/kernel_interface/firedrake_loopy.py b/tsfc/kernel_interface/firedrake_loopy.py deleted file mode 100644 index c6d671b7..00000000 --- a/tsfc/kernel_interface/firedrake_loopy.py +++ /dev/null @@ -1,435 +0,0 @@ -import numpy -from collections import namedtuple, OrderedDict -from functools import partial - -from ufl import Coefficient, FunctionSpace -from ufl.domain import extract_unique_domain -from finat.ufl import MixedElement as ufl_MixedElement, FiniteElement - -import gem -from gem.flop_count import count_flops - -import loopy as lp - -from tsfc import kernel_args -from tsfc.finatinterface import create_element -from tsfc.kernel_interface.common import KernelBuilderBase as _KernelBuilderBase, KernelBuilderMixin, get_index_names, check_requirements, prepare_coefficient, prepare_arguments, prepare_constant -from tsfc.loopy import generate as generate_loopy - - -# Expression kernel description type -ExpressionKernel = namedtuple('ExpressionKernel', ['ast', 'oriented', 'needs_cell_sizes', - 'coefficient_numbers', - 'needs_external_coords', - 'tabulations', 'name', 'arguments', - 'flop_count', 'event']) - - -def make_builder(*args, **kwargs): - return partial(KernelBuilder, *args, **kwargs) - - -class Kernel: - __slots__ = ("ast", "arguments", "integral_type", "oriented", "subdomain_id", - "domain_number", "needs_cell_sizes", "tabulations", - "coefficient_numbers", "name", "flop_count", "event", - "__weakref__") - """A compiled Kernel object. - - :kwarg ast: The loopy kernel object. - :kwarg integral_type: The type of integral. - :kwarg oriented: Does the kernel require cell_orientations. - :kwarg subdomain_id: What is the subdomain id for this kernel. - :kwarg domain_number: Which domain number in the original form - does this kernel correspond to (can be used to index into - original_form.ufl_domains() to get the correct domain). - :kwarg coefficient_numbers: A list of which coefficients from the - form the kernel needs. - :kwarg tabulations: The runtime tabulations this kernel requires - :kwarg needs_cell_sizes: Does the kernel require cell sizes. - :kwarg name: The name of this kernel. - :kwarg flop_count: Estimated total flops for this kernel. - :kwarg event: name for logging event - """ - def __init__(self, ast=None, arguments=None, integral_type=None, oriented=False, - subdomain_id=None, domain_number=None, - coefficient_numbers=(), - needs_cell_sizes=False, - tabulations=None, - flop_count=0, - name=None, - event=None): - # Defaults - self.ast = ast - self.arguments = arguments - self.integral_type = integral_type - self.oriented = oriented - self.domain_number = domain_number - self.subdomain_id = subdomain_id - self.coefficient_numbers = coefficient_numbers - self.needs_cell_sizes = needs_cell_sizes - self.tabulations = tabulations - self.flop_count = flop_count - self.name = name - self.event = event - - -class KernelBuilderBase(_KernelBuilderBase): - - def __init__(self, scalar_type, interior_facet=False): - """Initialise a kernel builder. - - :arg interior_facet: kernel accesses two cells - """ - super().__init__(scalar_type=scalar_type, interior_facet=interior_facet) - - # Cell orientation - if self.interior_facet: - cell_orientations = gem.Variable("cell_orientations", (2,)) - self._cell_orientations = (gem.Indexed(cell_orientations, (0,)), - gem.Indexed(cell_orientations, (1,))) - else: - cell_orientations = gem.Variable("cell_orientations", (1,)) - self._cell_orientations = (gem.Indexed(cell_orientations, (0,)),) - - def _coefficient(self, coefficient, name): - """Prepare a coefficient. Adds glue code for the coefficient - and adds the coefficient to the coefficient map. - - :arg coefficient: :class:`ufl.Coefficient` - :arg name: coefficient name - :returns: GEM expression representing the coefficient - """ - expr = prepare_coefficient(coefficient, name, interior_facet=self.interior_facet) - self.coefficient_map[coefficient] = expr - return expr - - def set_cell_sizes(self, domain): - """Setup a fake coefficient for "cell sizes". - - :arg domain: The domain of the integral. - - This is required for scaling of derivative basis functions on - physically mapped elements (Argyris, Bell, etc...). We need a - measure of the mesh size around each vertex (hence this lives - in P1). - - Should the domain have topological dimension 0 this does - nothing. - """ - if domain.ufl_cell().topological_dimension() > 0: - # Can't create P1 since only P0 is a valid finite element if - # topological_dimension is 0 and the concept of "cell size" - # is not useful for a vertex. - f = Coefficient(FunctionSpace(domain, FiniteElement("P", domain.ufl_cell(), 1))) - expr = prepare_coefficient(f, "cell_sizes", interior_facet=self.interior_facet) - self._cell_sizes = expr - - def create_element(self, element, **kwargs): - """Create a FInAT element (suitable for tabulating with) given - a UFL element.""" - return create_element(element, **kwargs) - - def generate_arg_from_variable(self, var, dtype=None): - """Generate kernel arg from a :class:`gem.Variable`. - - :arg var: a :class:`gem.Variable` - :arg dtype: dtype of the kernel arg - :returns: kernel arg - """ - return lp.GlobalArg(var.name, dtype=dtype or self.scalar_type, shape=var.shape) - - def generate_arg_from_expression(self, expr, dtype=None): - """Generate kernel arg from gem expression(s). - - :arg expr: gem expression(s) representing a coefficient or the output tensor - :arg dtype: dtype of the kernel arg - :returns: kernel arg - """ - var, = gem.extract_type(expr if isinstance(expr, tuple) else (expr, ), gem.Variable) - return self.generate_arg_from_variable(var, dtype=dtype or self.scalar_type) - - -class ExpressionKernelBuilder(KernelBuilderBase): - """Builds expression kernels for UFL interpolation in Firedrake.""" - - def __init__(self, scalar_type): - super(ExpressionKernelBuilder, self).__init__(scalar_type=scalar_type) - self.oriented = False - self.cell_sizes = False - - def set_coefficients(self, coefficients): - """Prepare the coefficients of the expression. - - :arg coefficients: UFL coefficients from Firedrake - """ - self.coefficient_split = {} - - for i, coefficient in enumerate(coefficients): - if type(coefficient.ufl_element()) == ufl_MixedElement: - subcoeffs = coefficient.subfunctions # Firedrake-specific - self.coefficient_split[coefficient] = subcoeffs - for j, subcoeff in enumerate(subcoeffs): - self._coefficient(subcoeff, f"w_{i}_{j}") - else: - self._coefficient(coefficient, f"w_{i}") - - def set_constants(self, constants): - for i, const in enumerate(constants): - gemexpr = prepare_constant(const, i) - self.constant_map[const] = gemexpr - - def set_coefficient_numbers(self, coefficient_numbers): - """Store the coefficient indices of the original form. - - :arg coefficient_numbers: Iterable of indices describing which coefficients - from the input expression need to be passed in to the kernel. - """ - self.coefficient_numbers = coefficient_numbers - - def register_requirements(self, ir): - """Inspect what is referenced by the IR that needs to be - provided by the kernel interface.""" - self.oriented, self.cell_sizes, self.tabulations = check_requirements(ir) - - def set_output(self, o): - """Produce the kernel return argument""" - loopy_arg = lp.GlobalArg(o.name, dtype=self.scalar_type, shape=o.shape) - self.output_arg = kernel_args.OutputKernelArg(loopy_arg) - - def construct_kernel(self, impero_c, index_names, needs_external_coords, log=False): - """Constructs an :class:`ExpressionKernel`. - - :arg impero_c: gem.ImperoC object that represents the kernel - :arg index_names: pre-assigned index names - :arg needs_external_coords: If ``True``, the first argument to - the kernel is an externally provided coordinate field. - :arg log: bool if the Kernel should be profiled with Log events - - :returns: :class:`ExpressionKernel` object - """ - args = [self.output_arg] - if self.oriented: - funarg = self.generate_arg_from_expression(self._cell_orientations, dtype=numpy.int32) - args.append(kernel_args.CellOrientationsKernelArg(funarg)) - if self.cell_sizes: - funarg = self.generate_arg_from_expression(self._cell_sizes) - args.append(kernel_args.CellSizesKernelArg(funarg)) - for _, expr in self.coefficient_map.items(): - # coefficient_map is OrderedDict. - funarg = self.generate_arg_from_expression(expr) - args.append(kernel_args.CoefficientKernelArg(funarg)) - - # now constants - for gemexpr in self.constant_map.values(): - funarg = self.generate_arg_from_expression(gemexpr) - args.append(kernel_args.ConstantKernelArg(funarg)) - - for name_, shape in self.tabulations: - tab_loopy_arg = lp.GlobalArg(name_, dtype=self.scalar_type, shape=shape) - args.append(kernel_args.TabulationKernelArg(tab_loopy_arg)) - - loopy_args = [arg.loopy_arg for arg in args] - - name = "expression_kernel" - loopy_kernel, event = generate_loopy(impero_c, loopy_args, self.scalar_type, - name, index_names, log=log) - return ExpressionKernel(loopy_kernel, self.oriented, self.cell_sizes, - self.coefficient_numbers, needs_external_coords, - self.tabulations, name, args, count_flops(impero_c), event) - - -class KernelBuilder(KernelBuilderBase, KernelBuilderMixin): - """Helper class for building a :class:`Kernel` object.""" - - def __init__(self, integral_data_info, scalar_type, - dont_split=(), diagonal=False): - """Initialise a kernel builder.""" - integral_type = integral_data_info.integral_type - super(KernelBuilder, self).__init__(scalar_type, integral_type.startswith("interior_facet")) - self.fem_scalar_type = scalar_type - - self.diagonal = diagonal - self.local_tensor = None - self.coefficient_split = {} - self.coefficient_number_index_map = OrderedDict() - self.dont_split = frozenset(dont_split) - - # Facet number - if integral_type in ['exterior_facet', 'exterior_facet_vert']: - facet = gem.Variable('facet', (1,)) - self._entity_number = {None: gem.VariableIndex(gem.Indexed(facet, (0,)))} - elif integral_type in ['interior_facet', 'interior_facet_vert']: - facet = gem.Variable('facet', (2,)) - self._entity_number = { - '+': gem.VariableIndex(gem.Indexed(facet, (0,))), - '-': gem.VariableIndex(gem.Indexed(facet, (1,))) - } - elif integral_type == 'interior_facet_horiz': - self._entity_number = {'+': 1, '-': 0} - - self.set_arguments(integral_data_info.arguments) - self.integral_data_info = integral_data_info - - def set_arguments(self, arguments): - """Process arguments. - - :arg arguments: :class:`ufl.Argument`s - :returns: GEM expression representing the return variable - """ - argument_multiindices = tuple(create_element(arg.ufl_element()).get_indices() - for arg in arguments) - if self.diagonal: - # Error checking occurs in the builder constructor. - # Diagonal assembly is obtained by using the test indices for - # the trial space as well. - a, _ = argument_multiindices - argument_multiindices = (a, a) - return_variables = prepare_arguments(arguments, - argument_multiindices, - interior_facet=self.interior_facet, - diagonal=self.diagonal) - self.return_variables = return_variables - self.argument_multiindices = argument_multiindices - - def set_coordinates(self, domain): - """Prepare the coordinate field. - - :arg domain: :class:`ufl.Domain` - """ - # Create a fake coordinate coefficient for a domain. - f = Coefficient(FunctionSpace(domain, domain.ufl_coordinate_element())) - self.domain_coordinate[domain] = f - self._coefficient(f, "coords") - - def set_coefficients(self, integral_data, form_data): - """Prepare the coefficients of the form. - - :arg integral_data: UFL integral data - :arg form_data: UFL form data - """ - # enabled_coefficients is a boolean array that indicates which - # of reduced_coefficients the integral requires. - n, k = 0, 0 - for i in range(len(integral_data.enabled_coefficients)): - if integral_data.enabled_coefficients[i]: - original = form_data.reduced_coefficients[i] - coefficient = form_data.function_replace_map[original] - if type(coefficient.ufl_element()) == ufl_MixedElement: - if original in self.dont_split: - self.coefficient_split[coefficient] = [coefficient] - self._coefficient(coefficient, f"w_{k}") - self.coefficient_number_index_map[coefficient] = (n, 0) - k += 1 - else: - self.coefficient_split[coefficient] = [] - for j, element in enumerate(coefficient.ufl_element().sub_elements): - c = Coefficient(FunctionSpace(extract_unique_domain(coefficient), element)) - self.coefficient_split[coefficient].append(c) - self._coefficient(c, f"w_{k}") - self.coefficient_number_index_map[c] = (n, j) - k += 1 - else: - self._coefficient(coefficient, f"w_{k}") - self.coefficient_number_index_map[coefficient] = (n, 0) - k += 1 - n += 1 - - def set_constants(self, constants): - for i, const in enumerate(constants): - gemexpr = prepare_constant(const, i) - self.constant_map[const] = gemexpr - - def register_requirements(self, ir): - """Inspect what is referenced by the IR that needs to be - provided by the kernel interface.""" - return check_requirements(ir) - - def construct_kernel(self, name, ctx, log=False): - """Construct a fully built :class:`Kernel`. - - This function contains the logic for building the argument - list for assembly kernels. - - :arg name: kernel name - :arg ctx: kernel builder context to get impero_c from - :arg log: bool if the Kernel should be profiled with Log events - :returns: :class:`Kernel` object - """ - impero_c, oriented, needs_cell_sizes, tabulations, active_variables = self.compile_gem(ctx) - if impero_c is None: - return self.construct_empty_kernel(name) - info = self.integral_data_info - # In the following funargs are only generated - # for gem expressions that are actually used; - # see `generate_arg_from_expression()` method. - # Specifically, funargs are not generated for - # unused components of mixed coefficients. - # Problem solving environment, such as Firedrake, - # will know which components have been included - # in the list of kernel arguments by investigating - # `Kernel.coefficient_numbers`. - # Add return arg - funarg = self.generate_arg_from_expression(self.return_variables) - args = [kernel_args.OutputKernelArg(funarg)] - # Add coordinates arg - coord = self.domain_coordinate[info.domain] - expr = self.coefficient_map[coord] - funarg = self.generate_arg_from_expression(expr) - args.append(kernel_args.CoordinatesKernelArg(funarg)) - if oriented: - funarg = self.generate_arg_from_expression(self._cell_orientations, dtype=numpy.int32) - args.append(kernel_args.CellOrientationsKernelArg(funarg)) - if needs_cell_sizes: - funarg = self.generate_arg_from_expression(self._cell_sizes) - args.append(kernel_args.CellSizesKernelArg(funarg)) - coefficient_indices = OrderedDict() - for coeff, (number, index) in self.coefficient_number_index_map.items(): - a = coefficient_indices.setdefault(number, []) - expr = self.coefficient_map[coeff] - var, = gem.extract_type(expr if isinstance(expr, tuple) else (expr, ), gem.Variable) - if var in active_variables: - funarg = self.generate_arg_from_expression(expr) - args.append(kernel_args.CoefficientKernelArg(funarg)) - a.append(index) - - # now constants - for gemexpr in self.constant_map.values(): - funarg = self.generate_arg_from_expression(gemexpr) - args.append(kernel_args.ConstantKernelArg(funarg)) - - coefficient_indices = tuple(tuple(v) for v in coefficient_indices.values()) - assert len(coefficient_indices) == len(info.coefficient_numbers) - if info.integral_type in ["exterior_facet", "exterior_facet_vert"]: - ext_loopy_arg = lp.GlobalArg("facet", numpy.uint32, shape=(1,)) - args.append(kernel_args.ExteriorFacetKernelArg(ext_loopy_arg)) - elif info.integral_type in ["interior_facet", "interior_facet_vert"]: - int_loopy_arg = lp.GlobalArg("facet", numpy.uint32, shape=(2,)) - args.append(kernel_args.InteriorFacetKernelArg(int_loopy_arg)) - for name_, shape in tabulations: - tab_loopy_arg = lp.GlobalArg(name_, dtype=self.scalar_type, shape=shape) - args.append(kernel_args.TabulationKernelArg(tab_loopy_arg)) - index_names = get_index_names(ctx['quadrature_indices'], self.argument_multiindices, ctx['index_cache']) - ast, event_name = generate_loopy(impero_c, [arg.loopy_arg for arg in args], - self.scalar_type, name, index_names, log=log) - flop_count = count_flops(impero_c) # Estimated total flops for this kernel. - return Kernel(ast=ast, - arguments=tuple(args), - integral_type=info.integral_type, - subdomain_id=info.subdomain_id, - domain_number=info.domain_number, - coefficient_numbers=tuple(zip(info.coefficient_numbers, coefficient_indices)), - oriented=oriented, - needs_cell_sizes=needs_cell_sizes, - tabulations=tabulations, - flop_count=flop_count, - name=name, - event=event_name) - - def construct_empty_kernel(self, name): - """Return None, since Firedrake needs no empty kernels. - - :arg name: function name - :returns: None - """ - return None diff --git a/tsfc/logging.py b/tsfc/logging.py deleted file mode 100644 index 6bced525..00000000 --- a/tsfc/logging.py +++ /dev/null @@ -1,6 +0,0 @@ -"""Logging for TSFC.""" - -import logging - -logger = logging.getLogger('tsfc') -logger.addHandler(logging.StreamHandler()) diff --git a/tsfc/loopy.py b/tsfc/loopy.py deleted file mode 100644 index 943143ff..00000000 --- a/tsfc/loopy.py +++ /dev/null @@ -1,549 +0,0 @@ -"""Generate loopy kernel from ImperoC tuple data. - -This is the final stage of code generation in TSFC.""" - -import numpy -from functools import singledispatch -from collections import defaultdict, OrderedDict - -from gem import gem, impero as imp -from gem.node import Memoizer - -import islpy as isl -import loopy as lp - -import pymbolic.primitives as p -from loopy.symbolic import SubArrayRef - -from pytools import UniqueNameGenerator - -from tsfc.parameters import is_complex - -from contextlib import contextmanager -from tsfc.parameters import target - - -def profile_insns(kernel_name, instructions, log=False): - if log: - event_name = "Log_Event_" + kernel_name - event_id_var_name = "ID_" + event_name - # Logging registration - # The events are registered in PyOP2 and the event id is passed onto the dll - preamble = "PetscLogEvent "+event_id_var_name+" = -1;" - # Profiling - prepend = [lp.CInstruction("", "PetscLogEventBegin("+event_id_var_name+",0,0,0,0);")] - append = [lp.CInstruction("", "PetscLogEventEnd("+event_id_var_name+",0,0,0,0);")] - instructions = prepend + instructions + append - return instructions, event_name, [(str(2**31-1)+"_"+kernel_name, preamble)] - else: - return instructions, None, None - - -@singledispatch -def _assign_dtype(expression, self): - return set.union(*map(self, expression.children)) - - -@_assign_dtype.register(gem.Terminal) -def _assign_dtype_terminal(expression, self): - return {self.scalar_type} - - -@_assign_dtype.register(gem.Zero) -@_assign_dtype.register(gem.Identity) -@_assign_dtype.register(gem.Delta) -def _assign_dtype_real(expression, self): - return {self.real_type} - - -@_assign_dtype.register(gem.Literal) -def _assign_dtype_identity(expression, self): - return {expression.array.dtype} - - -@_assign_dtype.register(gem.Power) -def _assign_dtype_power(expression, self): - # Conservative - return {self.scalar_type} - - -@_assign_dtype.register(gem.MathFunction) -def _assign_dtype_mathfunction(expression, self): - if expression.name in {"abs", "real", "imag"}: - return {self.real_type} - elif expression.name == "sqrt": - return {self.scalar_type} - else: - return set.union(*map(self, expression.children)) - - -@_assign_dtype.register(gem.MinValue) -@_assign_dtype.register(gem.MaxValue) -def _assign_dtype_minmax(expression, self): - # UFL did correctness checking - return {self.real_type} - - -@_assign_dtype.register(gem.Conditional) -def _assign_dtype_conditional(expression, self): - return set.union(*map(self, expression.children[1:])) - - -@_assign_dtype.register(gem.Comparison) -@_assign_dtype.register(gem.LogicalNot) -@_assign_dtype.register(gem.LogicalAnd) -@_assign_dtype.register(gem.LogicalOr) -def _assign_dtype_logical(expression, self): - return {numpy.int8} - - -def assign_dtypes(expressions, scalar_type): - """Assign numpy data types to expressions. - - Used for declaring temporaries when converting from Impero to lower level code. - - :arg expressions: List of GEM expressions. - :arg scalar_type: Default scalar type. - - :returns: list of tuples (expression, dtype).""" - mapper = Memoizer(_assign_dtype) - mapper.scalar_type = scalar_type - mapper.real_type = numpy.finfo(scalar_type).dtype - return [(e, numpy.result_type(*mapper(e))) for e in expressions] - - -class LoopyContext(object): - def __init__(self, target=None): - self.indices = {} # indices for declarations and referencing values, from ImperoC - self.active_indices = {} # gem index -> pymbolic variable - self.index_extent = OrderedDict() # pymbolic variable for indices -> extent - self.gem_to_pymbolic = {} # gem node -> pymbolic variable - self.name_gen = UniqueNameGenerator() - self.target = target - - def fetch_multiindex(self, multiindex): - indices = [] - for index in multiindex: - if isinstance(index, gem.Index): - indices.append(self.active_indices[index]) - elif isinstance(index, gem.VariableIndex): - indices.append(expression(index.expression, self)) - else: - assert isinstance(index, int) - indices.append(index) - return tuple(indices) - - # Generate index from gem multiindex - def gem_to_pym_multiindex(self, multiindex): - indices = [] - for index in multiindex: - assert index.extent - if not index.name: - name = self.name_gen(self.index_names[index]) - else: - name = index.name - self.index_extent[name] = index.extent - indices.append(p.Variable(name)) - return tuple(indices) - - # Generate index from shape - def pymbolic_multiindex(self, shape): - indices = [] - for extent in shape: - name = self.name_gen(self.index_names[extent]) - self.index_extent[name] = extent - indices.append(p.Variable(name)) - return tuple(indices) - - # Generate pym variable from gem - def pymbolic_variable_and_destruct(self, node): - pym = self.pymbolic_variable(node) - if isinstance(pym, p.Subscript): - return pym.aggregate, pym.index_tuple - else: - return pym, () - - # Generate pym variable or subscript - def pymbolic_variable(self, node): - pym = self._gem_to_pym_var(node) - if node in self.indices: - indices = self.fetch_multiindex(self.indices[node]) - if indices: - return p.Subscript(pym, indices) - return pym - - def _gem_to_pym_var(self, node): - try: - pym = self.gem_to_pymbolic[node] - except KeyError: - name = self.name_gen(node.name) - pym = p.Variable(name) - self.gem_to_pymbolic[node] = pym - return pym - - def active_inames(self): - # Return all active indices - return frozenset([i.name for i in self.active_indices.values()]) - - -@contextmanager -def active_indices(mapping, ctx): - """Push active indices onto context. - :arg mapping: dict mapping gem indices to pymbolic index expressions - :arg ctx: code generation context. - :returns: new code generation context.""" - ctx.active_indices.update(mapping) - yield ctx - for key in mapping: - ctx.active_indices.pop(key) - - -def generate(impero_c, args, scalar_type, kernel_name="loopy_kernel", index_names=[], - return_increments=True, log=False): - """Generates loopy code. - - :arg impero_c: ImperoC tuple with Impero AST and other data - :arg args: list of loopy.GlobalArgs - :arg scalar_type: type of scalars as C typename string - :arg kernel_name: function name of the kernel - :arg index_names: pre-assigned index names - :arg return_increments: Does codegen for Return nodes increment the lvalue, or assign? - :arg log: bool if the Kernel should be profiled with Log events - :returns: loopy kernel - """ - ctx = LoopyContext(target=target) - ctx.indices = impero_c.indices - ctx.index_names = defaultdict(lambda: "i", index_names) - ctx.epsilon = numpy.finfo(scalar_type).resolution - ctx.scalar_type = scalar_type - ctx.return_increments = return_increments - - # Create arguments - data = list(args) - for i, (temp, dtype) in enumerate(assign_dtypes(impero_c.temporaries, scalar_type)): - name = "t%d" % i - if isinstance(temp, gem.Constant): - data.append(lp.TemporaryVariable(name, shape=temp.shape, dtype=dtype, initializer=temp.array, address_space=lp.AddressSpace.LOCAL, read_only=True)) - else: - shape = tuple([i.extent for i in ctx.indices[temp]]) + temp.shape - data.append(lp.TemporaryVariable(name, shape=shape, dtype=dtype, initializer=None, address_space=lp.AddressSpace.LOCAL, read_only=False)) - ctx.gem_to_pymbolic[temp] = p.Variable(name) - - # Create instructions - instructions = statement(impero_c.tree, ctx) - - # add a no-op touching all kernel arguments to make sure they - # are not silently dropped - noop = lp.CInstruction( - (), "", read_variables=frozenset({a.name for a in args}), - within_inames=frozenset(), within_inames_is_final=True) - instructions.append(noop) - - # Profile the instructions - instructions, event_name, preamble = profile_insns(kernel_name, instructions, log) - - # Create domains - domains = create_domains(ctx.index_extent.items()) - - # Create loopy kernel - knl = lp.make_kernel( - domains, - instructions, - data, - name=kernel_name, - target=target, - seq_dependencies=True, - silenced_warnings=["summing_if_branches_ops"], - lang_version=(2018, 2), - preambles=preamble - ) - - # Prevent loopy interchange by loopy - knl = lp.prioritize_loops(knl, ",".join(ctx.index_extent.keys())) - - return knl, event_name - - -def create_domains(indices): - """ Create ISL domains from indices - - :arg indices: iterable of (index_name, extent) pairs - :returns: A list of ISL sets representing the iteration domain of the indices.""" - - domains = [] - for idx, extent in indices: - inames = isl.make_zero_and_vars([idx]) - domains.append(((inames[0].le_set(inames[idx])) & (inames[idx].lt_set(inames[0] + extent)))) - - if not domains: - domains = [isl.BasicSet("[] -> {[]}")] - return domains - - -@singledispatch -def statement(tree, ctx): - """Translates an Impero (sub)tree into a loopy instructions corresponding - to a C statement. - - :arg tree: Impero (sub)tree - :arg ctx: miscellaneous code generation data - :returns: list of loopy instructions - """ - raise AssertionError("cannot generate loopy from %s" % type(tree)) - - -@statement.register(imp.Block) -def statement_block(tree, ctx): - from itertools import chain - return list(chain(*(statement(child, ctx) for child in tree.children))) - - -@statement.register(imp.For) -def statement_for(tree, ctx): - extent = tree.index.extent - assert extent - idx = ctx.name_gen(ctx.index_names[tree.index]) - ctx.index_extent[idx] = extent - with active_indices({tree.index: p.Variable(idx)}, ctx) as ctx_active: - return statement(tree.children[0], ctx_active) - - -@statement.register(imp.Initialise) -def statement_initialise(leaf, ctx): - return [lp.Assignment(expression(leaf.indexsum, ctx), 0.0, within_inames=ctx.active_inames())] - - -@statement.register(imp.Accumulate) -def statement_accumulate(leaf, ctx): - lhs = expression(leaf.indexsum, ctx) - rhs = lhs + expression(leaf.indexsum.children[0], ctx) - return [lp.Assignment(lhs, rhs, within_inames=ctx.active_inames())] - - -@statement.register(imp.Return) -def statement_return(leaf, ctx): - lhs = expression(leaf.variable, ctx) - rhs = expression(leaf.expression, ctx) - if ctx.return_increments: - rhs = lhs + rhs - return [lp.Assignment(lhs, rhs, within_inames=ctx.active_inames())] - - -@statement.register(imp.ReturnAccumulate) -def statement_returnaccumulate(leaf, ctx): - lhs = expression(leaf.variable, ctx) - rhs = lhs + expression(leaf.indexsum.children[0], ctx) - return [lp.Assignment(lhs, rhs, within_inames=ctx.active_inames())] - - -@statement.register(imp.Evaluate) -def statement_evaluate(leaf, ctx): - expr = leaf.expression - if isinstance(expr, gem.ListTensor): - ops = [] - var, index = ctx.pymbolic_variable_and_destruct(expr) - for multiindex, value in numpy.ndenumerate(expr.array): - ops.append(lp.Assignment(p.Subscript(var, index + multiindex), expression(value, ctx), within_inames=ctx.active_inames())) - return ops - elif isinstance(expr, gem.Constant): - return [] - elif isinstance(expr, gem.ComponentTensor): - idx = ctx.gem_to_pym_multiindex(expr.multiindex) - var, sub_idx = ctx.pymbolic_variable_and_destruct(expr) - lhs = p.Subscript(var, idx + sub_idx) - with active_indices(dict(zip(expr.multiindex, idx)), ctx) as ctx_active: - return [lp.Assignment(lhs, expression(expr.children[0], ctx_active), within_inames=ctx_active.active_inames())] - elif isinstance(expr, gem.Inverse): - idx = ctx.pymbolic_multiindex(expr.shape) - var = ctx.pymbolic_variable(expr) - lhs = (SubArrayRef(idx, p.Subscript(var, idx)),) - - idx_reads = ctx.pymbolic_multiindex(expr.children[0].shape) - var_reads = ctx.pymbolic_variable(expr.children[0]) - reads = (SubArrayRef(idx_reads, p.Subscript(var_reads, idx_reads)),) - rhs = p.Call(p.Variable("inverse"), reads) - - return [lp.CallInstruction(lhs, rhs, within_inames=ctx.active_inames())] - elif isinstance(expr, gem.Solve): - idx = ctx.pymbolic_multiindex(expr.shape) - var = ctx.pymbolic_variable(expr) - lhs = (SubArrayRef(idx, p.Subscript(var, idx)),) - - reads = [] - for child in expr.children: - idx_reads = ctx.pymbolic_multiindex(child.shape) - var_reads = ctx.pymbolic_variable(child) - reads.append(SubArrayRef(idx_reads, p.Subscript(var_reads, idx_reads))) - rhs = p.Call(p.Variable("solve"), tuple(reads)) - - return [lp.CallInstruction(lhs, rhs, within_inames=ctx.active_inames())] - else: - return [lp.Assignment(ctx.pymbolic_variable(expr), expression(expr, ctx, top=True), within_inames=ctx.active_inames())] - - -def expression(expr, ctx, top=False): - """Translates GEM expression into a pymbolic expression - - :arg expr: GEM expression - :arg ctx: miscellaneous code generation data - :arg top: do not generate temporary reference for the root node - :returns: pymbolic expression - """ - if not top and expr in ctx.gem_to_pymbolic: - return ctx.pymbolic_variable(expr) - else: - return _expression(expr, ctx) - - -@singledispatch -def _expression(expr, ctx): - raise AssertionError("cannot generate expression from %s" % type(expr)) - - -@_expression.register(gem.Failure) -def _expression_failure(expr, ctx): - raise expr.exception - - -@_expression.register(gem.Product) -def _expression_product(expr, ctx): - return p.Product(tuple(expression(c, ctx) for c in expr.children)) - - -@_expression.register(gem.Sum) -def _expression_sum(expr, ctx): - return p.Sum(tuple(expression(c, ctx) for c in expr.children)) - - -@_expression.register(gem.Division) -def _expression_division(expr, ctx): - return p.Quotient(*(expression(c, ctx) for c in expr.children)) - - -@_expression.register(gem.Power) -def _expression_power(expr, ctx): - return p.Variable("pow")(*(expression(c, ctx) for c in expr.children)) - - -@_expression.register(gem.MathFunction) -def _expression_mathfunction(expr, ctx): - if expr.name.startswith('cyl_bessel_'): - # Bessel functions - if is_complex(ctx.scalar_type): - raise NotImplementedError("Bessel functions for complex numbers: " - "missing implementation") - nu, arg = expr.children - nu_ = expression(nu, ctx) - arg_ = expression(arg, ctx) - if isinstance(ctx.target, lp.target.c.CWithGNULibcTarget): - # Generate right functions calls to gnulibc bessel functions - # cyl_bessel_{jy} -> bessel_{jy} - name = expr.name[4:] - return p.Variable(f"{name}n")(int(nu_), arg_) - else: - # Modified Bessel functions (C++ only) - # These mappings work for FEniCS only, and fail with Firedrake - # since no Boost available. - # Is this actually still supported/has ever been used by anyone? - if expr.name in {'cyl_bessel_i', 'cyl_bessel_k'}: - name = 'boost::math::' + expr.name - return p.Variable(name)(nu_, arg_) - else: - # cyl_bessel_{jy} -> {jy} - name = expr.name[-1:] - if nu == gem.Zero(): - return p.Variable(f"{name}0")(arg_) - elif nu == gem.one: - return p.Variable(f"{name}1")(arg_) - else: - return p.Variable(f"{name}n")(nu_, arg_) - else: - if expr.name == "ln": - name = "log" - else: - name = expr.name - # Not all mathfunctions apply to complex numbers, but this - # will be picked up in loopy. This way we allow erf(real(...)) - # in complex mode (say). - return p.Variable(name)(*(expression(c, ctx) for c in expr.children)) - - -@_expression.register(gem.MinValue) -def _expression_minvalue(expr, ctx): - return p.Variable("min")(*[expression(c, ctx) for c in expr.children]) - - -@_expression.register(gem.MaxValue) -def _expression_maxvalue(expr, ctx): - return p.Variable("max")(*[expression(c, ctx) for c in expr.children]) - - -@_expression.register(gem.Comparison) -def _expression_comparison(expr, ctx): - left, right = [expression(c, ctx) for c in expr.children] - return p.Comparison(left, expr.operator, right) - - -@_expression.register(gem.LogicalNot) -def _expression_logicalnot(expr, ctx): - child, = expr.children - return p.LogicalNot(expression(child, ctx)) - - -@_expression.register(gem.LogicalAnd) -def _expression_logicaland(expr, ctx): - return p.LogicalAnd(tuple([expression(c, ctx) for c in expr.children])) - - -@_expression.register(gem.LogicalOr) -def _expression_logicalor(expr, ctx): - return p.LogicalOr(tuple([expression(c, ctx) for c in expr.children])) - - -@_expression.register(gem.Conditional) -def _expression_conditional(expr, ctx): - return p.If(*[expression(c, ctx) for c in expr.children]) - - -@_expression.register(gem.Constant) -def _expression_scalar(expr, ctx): - assert not expr.shape - v = expr.value - if numpy.isnan(v): - return p.Variable("NAN") - r = numpy.round(v, 1) - if r and numpy.abs(v - r) < ctx.epsilon: - return r - return v - - -@_expression.register(gem.Variable) -def _expression_variable(expr, ctx): - return ctx.pymbolic_variable(expr) - - -@_expression.register(gem.Indexed) -def _expression_indexed(expr, ctx): - rank = ctx.fetch_multiindex(expr.multiindex) - var = expression(expr.children[0], ctx) - if isinstance(var, p.Subscript): - rank = var.index + rank - var = var.aggregate - return p.Subscript(var, rank) - - -@_expression.register(gem.FlexiblyIndexed) -def _expression_flexiblyindexed(expr, ctx): - var = expression(expr.children[0], ctx) - - rank = [] - for off, idxs in expr.dim2idxs: - for index, stride in idxs: - assert isinstance(index, gem.Index) - - rank_ = [off] - for index, stride in idxs: - rank_.append(p.Product((ctx.active_indices[index], stride))) - rank.append(p.Sum(tuple(rank_))) - - return p.Subscript(var, tuple(rank)) diff --git a/tsfc/modified_terminals.py b/tsfc/modified_terminals.py deleted file mode 100644 index 8c5162bf..00000000 --- a/tsfc/modified_terminals.py +++ /dev/null @@ -1,179 +0,0 @@ -# -*- coding: utf-8 -*- -# Copyright (C) 2011-2015 Martin Sandve Alnæs -# -# This file is part of UFLACS. -# -# UFLACS 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. -# -# UFLACS 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 UFLACS. If not, see . -# -# Modified by Miklós Homolya, 2016. - -"""Definitions of 'modified terminals', a core concept in uflacs.""" - -from ufl.classes import (ReferenceValue, ReferenceGrad, - NegativeRestricted, PositiveRestricted, - Restricted, ConstantValue, - Jacobian, SpatialCoordinate, Zero) -from ufl.checks import is_cellwise_constant -from ufl.domain import extract_unique_domain - - -class ModifiedTerminal(object): - - """A modified terminal expression is an object of a Terminal subtype, wrapped in terminal modifier types. - - The variables of this class are: - - expr - The original UFL expression - - terminal - the underlying Terminal object - local_derivatives - tuple of ints, each meaning derivative in that local direction - reference_value - bool, whether this is represented in reference frame - restriction - None, '+' or '-' - """ - - def __init__(self, expr, terminal, local_derivatives, restriction, reference_value): - # The original expression - self.expr = expr - - # The underlying terminal expression - self.terminal = terminal - - # Components - self.reference_value = reference_value - self.restriction = restriction - - # Derivatives - self.local_derivatives = local_derivatives - - def as_tuple(self): - t = self.terminal - rv = self.reference_value - ld = self.local_derivatives - r = self.restriction - return (t, rv, ld, r) - - def __hash__(self): - return hash(self.as_tuple()) - - def __eq__(self, other): - return isinstance(other, ModifiedTerminal) and self.as_tuple() == other.as_tuple() - - def __lt__(self, other): - return self.as_tuple() < other.as_tuple() - - def __str__(self): - s = [] - s += ["terminal: {0}".format(self.terminal)] - s += ["local_derivatives: {0}".format(self.local_derivatives)] - s += ["restriction: {0}".format(self.restriction)] - return '\n'.join(s) - - -def is_modified_terminal(v): - "Check if v is a terminal or a terminal wrapped in terminal modifier types." - while not v._ufl_is_terminal_: - if v._ufl_is_terminal_modifier_: - v = v.ufl_operands[0] - else: - return False - return True - - -def strip_modified_terminal(v): - "Extract core Terminal from a modified terminal or return None." - while not v._ufl_is_terminal_: - if v._ufl_is_terminal_modifier_: - v = v.ufl_operands[0] - else: - return None - return v - - -def analyse_modified_terminal(expr): - """Analyse a so-called 'modified terminal' expression and return its properties in more compact form. - - A modified terminal expression is an object of a Terminal subtype, wrapped in terminal modifier types. - - The wrapper types can include 0-* Grad or ReferenceGrad objects, - and 0-1 ReferenceValue, 0-1 Restricted, and 0-1 Indexed objects. - """ - # Data to determine - local_derivatives = 0 - reference_value = None - restriction = None - - # Start with expr and strip away layers of modifiers - t = expr - while not t._ufl_is_terminal_: - if isinstance(t, ReferenceValue): - assert reference_value is None, "Got twice pulled back terminal!" - reference_value = True - t, = t.ufl_operands - - elif isinstance(t, ReferenceGrad): - local_derivatives += 1 - t, = t.ufl_operands - - elif isinstance(t, Restricted): - assert restriction is None, "Got twice restricted terminal!" - restriction = t._side - t, = t.ufl_operands - - elif t._ufl_terminal_modifiers_: - raise ValueError("Missing handler for terminal modifier type %s, object is %s." % (type(t), repr(t))) - - else: - raise ValueError("Unexpected type %s object %s." % (type(t), repr(t))) - - # Make reference_value true or false - if reference_value is None: - reference_value = False - - # Consistency check - if isinstance(t, (SpatialCoordinate, Jacobian)): - pass - else: - if local_derivatives and not reference_value: - raise ValueError("Local derivatives of non-local value?") - - return ModifiedTerminal(expr, t, local_derivatives, restriction, reference_value) - - -def construct_modified_terminal(mt, terminal): - """Construct a modified terminal given terminal modifiers from an - analysed modified terminal and a terminal.""" - expr = terminal - - if mt.reference_value: - expr = ReferenceValue(expr) - - dim = extract_unique_domain(expr).topological_dimension() - for n in range(mt.local_derivatives): - # Return zero if expression is trivially constant. This has to - # happen here because ReferenceGrad has no access to the - # topological dimension of a literal zero. - if is_cellwise_constant(expr): - expr = Zero(expr.ufl_shape + (dim,), expr.ufl_free_indices, - expr.ufl_index_dimensions) - else: - expr = ReferenceGrad(expr) - - # No need to apply restrictions to ConstantValue terminals - if not isinstance(expr, ConstantValue): - if mt.restriction == '+': - expr = PositiveRestricted(expr) - elif mt.restriction == '-': - expr = NegativeRestricted(expr) - - return expr diff --git a/tsfc/parameters.py b/tsfc/parameters.py deleted file mode 100644 index 1277713a..00000000 --- a/tsfc/parameters.py +++ /dev/null @@ -1,36 +0,0 @@ -import numpy -from loopy.target.c import CWithGNULibcTarget - - -PARAMETERS = { - "quadrature_rule": "auto", - "quadrature_degree": "auto", - - # Default mode - "mode": "spectral", - - # Maximum extent to unroll index sums. Default is 3, so that loops - # over geometric dimensions are unrolled; this improves assembly - # performance. Can be disabled by setting it to None, False or 0; - # that makes compilation time much shorter. - "unroll_indexsum": 3, - - # Scalar type numpy dtype - "scalar_type": numpy.dtype(numpy.float64), - - # So that tests pass (needs to match scalar_type) - "scalar_type_c": "double", -} - - -target = CWithGNULibcTarget() - - -def default_parameters(): - return PARAMETERS.copy() - - -def is_complex(scalar_type): - """Decides complex mode based on scalar type.""" - return scalar_type and (isinstance(scalar_type, numpy.dtype) and scalar_type.kind == 'c') \ - or (isinstance(scalar_type, str) and "complex" in scalar_type) diff --git a/tsfc/spectral.py b/tsfc/spectral.py deleted file mode 100644 index 8ee83848..00000000 --- a/tsfc/spectral.py +++ /dev/null @@ -1,193 +0,0 @@ -from collections import OrderedDict, defaultdict, namedtuple -from functools import partial, reduce -from itertools import chain, zip_longest - -from gem.gem import Delta, Indexed, Sum, index_sum, one -from gem.node import Memoizer -from gem.optimise import delta_elimination as _delta_elimination -from gem.optimise import remove_componenttensors, replace_division, unroll_indexsum -from gem.refactorise import ATOMIC, COMPOUND, OTHER, MonomialSum, collect_monomials -from gem.unconcatenate import unconcatenate -from gem.coffee import optimise_monomial_sum -from gem.utils import groupby - - -Integral = namedtuple('Integral', ['expression', - 'quadrature_multiindex', - 'argument_indices']) - - -def Integrals(expressions, quadrature_multiindex, argument_multiindices, parameters): - """Constructs an integral representation for each GEM integrand - expression. - - :arg expressions: integrand multiplied with quadrature weight; - multi-root GEM expression DAG - :arg quadrature_multiindex: quadrature multiindex (tuple) - :arg argument_multiindices: tuple of argument multiindices, - one multiindex for each argument - :arg parameters: parameters dictionary - - :returns: list of integral representations - """ - # Rewrite: a / b => a * (1 / b) - expressions = replace_division(expressions) - - # Unroll - max_extent = parameters["unroll_indexsum"] - if max_extent: - def predicate(index): - return index.extent <= max_extent - expressions = unroll_indexsum(expressions, predicate=predicate) - - expressions = [index_sum(e, quadrature_multiindex) for e in expressions] - argument_indices = tuple(chain(*argument_multiindices)) - return [Integral(e, quadrature_multiindex, argument_indices) for e in expressions] - - -def _delta_inside(node, self): - """Does node contain a Delta?""" - return any(isinstance(child, Delta) or self(child) - for child in node.children) - - -def flatten(var_reps, index_cache): - quadrature_indices = OrderedDict() - - pairs = [] # assignment pairs - for variable, reps in var_reps: - # Extract argument indices - argument_indices, = set(r.argument_indices for r in reps) - assert set(variable.free_indices) == set(argument_indices) - - # Extract and verify expressions - expressions = [r.expression for r in reps] - assert all(set(e.free_indices) <= set(argument_indices) - for e in expressions) - - # Save assignment pair - pairs.append((variable, reduce(Sum, expressions))) - - # Collect quadrature_indices - for r in reps: - quadrature_indices.update(zip_longest(r.quadrature_multiindex, ())) - - # Split Concatenate nodes - pairs = unconcatenate(pairs, cache=index_cache) - - def group_key(pair): - variable, expression = pair - return frozenset(variable.free_indices) - - delta_inside = Memoizer(_delta_inside) - # Variable ordering after delta cancellation - narrow_variables = OrderedDict() - # Assignments are variable -> MonomialSum map - delta_simplified = defaultdict(MonomialSum) - # Group assignment pairs by argument indices - for free_indices, pair_group in groupby(pairs, group_key): - variables, expressions = zip(*pair_group) - # Argument factorise expressions - classifier = partial(classify, set(free_indices), delta_inside=delta_inside) - monomial_sums = collect_monomials(expressions, classifier) - # For each monomial, apply delta cancellation and insert - # result into delta_simplified. - for variable, monomial_sum in zip(variables, monomial_sums): - for monomial in monomial_sum: - var, s, a, r = delta_elimination(variable, *monomial) - narrow_variables.setdefault(var) - delta_simplified[var].add(s, a, r) - - # Final factorisation - for variable in narrow_variables: - monomial_sum = delta_simplified[variable] - # Collect sum indices applicable to the current MonomialSum - sum_indices = set().union(*[m.sum_indices for m in monomial_sum]) - # Put them in a deterministic order - sum_indices = [i for i in quadrature_indices if i in sum_indices] - # Sort for increasing index extent, this obtains the good - # factorisation for triangle x interval cells. Python sort is - # stable, so in the common case when index extents are equal, - # the previous deterministic ordering applies which is good - # for getting smaller temporaries. - sum_indices = sorted(sum_indices, key=lambda index: index.extent) - # Apply sum factorisation combined with COFFEE technology - expression = sum_factorise(variable, sum_indices, monomial_sum) - yield (variable, expression) - - -finalise_options = dict(replace_delta=False) - - -def classify(argument_indices, expression, delta_inside): - """Classifier for argument factorisation""" - n = len(argument_indices.intersection(expression.free_indices)) - if n == 0: - return OTHER - elif n == 1: - if isinstance(expression, (Delta, Indexed)) and not delta_inside(expression): - return ATOMIC - else: - return COMPOUND - else: - return COMPOUND - - -def delta_elimination(variable, sum_indices, args, rest): - """IndexSum-Delta cancellation for monomials.""" - factors = list(args) + [variable, rest] # construct factors - - def prune(factors): - # Skip last factor (``rest``, see above) which can be - # arbitrarily complicated, so its pruning may be expensive, - # and its early pruning brings no advantages. - result = remove_componenttensors(factors[:-1]) - result.append(factors[-1]) - return result - - # Cancel sum indices - sum_indices, factors = _delta_elimination(sum_indices, factors) - factors = prune(factors) - - # Cancel variable indices - var_indices, factors = _delta_elimination(variable.free_indices, factors) - factors = prune(factors) - - # Destructure factors after cancellation - rest = factors.pop() - variable = factors.pop() - args = [f for f in factors if f != one] - - assert set(var_indices) == set(variable.free_indices) - return variable, sum_indices, args, rest - - -def sum_factorise(variable, tail_ordering, monomial_sum): - if tail_ordering: - key_ordering = OrderedDict() - sub_monosums = defaultdict(MonomialSum) - for sum_indices, atomics, rest in monomial_sum: - # Pull out those sum indices that are not contained in the - # tail ordering, together with those atomics which do not - # share free indices with the tail ordering. - # - # Based on this, split the monomial sum, then recursively - # optimise each sub monomial sum with the first tail index - # removed. - tail_indices = tuple(i for i in sum_indices if i in tail_ordering) - tail_atomics = tuple(a for a in atomics - if set(tail_indices) & set(a.free_indices)) - head_indices = tuple(i for i in sum_indices if i not in tail_ordering) - head_atomics = tuple(a for a in atomics if a not in tail_atomics) - key = (head_indices, head_atomics) - key_ordering.setdefault(key) - sub_monosums[key].add(tail_indices, tail_atomics, rest) - sub_monosums = [(k, sub_monosums[k]) for k in key_ordering] - - monomial_sum = MonomialSum() - for (sum_indices, atomics), monosum in sub_monosums: - new_rest = sum_factorise(variable, tail_ordering[1:], monosum) - monomial_sum.add(sum_indices, atomics, new_rest) - - # Use COFFEE algorithm to optimise the monomial sum - return optimise_monomial_sum(monomial_sum, variable.index_ordering()) diff --git a/tsfc/tensor.py b/tsfc/tensor.py deleted file mode 100644 index 8902b470..00000000 --- a/tsfc/tensor.py +++ /dev/null @@ -1,93 +0,0 @@ -from collections import defaultdict -from functools import partial, reduce -from itertools import count - -import numpy - -import gem -from gem.optimise import remove_componenttensors, unroll_indexsum -from gem.refactorise import ATOMIC, COMPOUND, OTHER, collect_monomials -from gem.unconcatenate import flatten as concatenate - - -def einsum(factors, sum_indices): - """Evaluates a tensor product at compile time. - - :arg factors: iterable of indexed GEM literals - :arg sum_indices: indices to sum over - :returns: a single indexed GEM literal - """ - # Maps the problem onto numpy.einsum - index2letter = defaultdict(partial(lambda c: chr(ord('i') + next(c)), count())) - operands = [] - subscript_parts = [] - for factor in factors: - literal, = factor.children - selectors = [] - letters = [] - for index in factor.multiindex: - if isinstance(index, int): - selectors.append(index) - else: - selectors.append(slice(None)) - letters.append(index2letter[index]) - operands.append(literal.array.__getitem__(tuple(selectors))) - subscript_parts.append(''.join(letters)) - - result_pairs = sorted((letter, index) - for index, letter in index2letter.items() - if index not in sum_indices) - - subscripts = ','.join(subscript_parts) + '->' + ''.join(l for l, i in result_pairs) - tensor = numpy.einsum(subscripts, *operands) - return gem.Indexed(gem.Literal(tensor), tuple(i for l, i in result_pairs)) - - -def Integrals(expressions, quadrature_multiindex, argument_multiindices, parameters): - # Concatenate - expressions = concatenate(expressions) - - # Unroll - max_extent = parameters["unroll_indexsum"] - if max_extent: - def predicate(index): - return index.extent <= max_extent - expressions = unroll_indexsum(expressions, predicate=predicate) - - # Refactorise - def classify(quadrature_indices, expression): - if not quadrature_indices.intersection(expression.free_indices): - return OTHER - elif isinstance(expression, gem.Indexed) and isinstance(expression.children[0], gem.Literal): - return ATOMIC - else: - return COMPOUND - classifier = partial(classify, set(quadrature_multiindex)) - - result = [] - for expr, monomial_sum in zip(expressions, collect_monomials(expressions, classifier)): - # Select quadrature indices that are present - quadrature_indices = set(index for index in quadrature_multiindex - if index in expr.free_indices) - - products = [] - for sum_indices, factors, rest in monomial_sum: - # Collapse quadrature literals for each monomial - if factors or quadrature_indices: - replacement = einsum(remove_componenttensors(factors), quadrature_indices) - else: - replacement = gem.Literal(1) - # Rebuild expression - products.append(gem.IndexSum(gem.Product(replacement, rest), sum_indices)) - result.append(reduce(gem.Sum, products, gem.Zero())) - return result - - -def flatten(var_reps, index_cache): - for variable, reps in var_reps: - expressions = reps - for expression in expressions: - yield (variable, expression) - - -finalise_options = {} diff --git a/tsfc/ufl2gem.py b/tsfc/ufl2gem.py deleted file mode 100644 index e283fb14..00000000 --- a/tsfc/ufl2gem.py +++ /dev/null @@ -1,160 +0,0 @@ -"""Translation of UFL tensor-algebra into GEM tensor-algebra.""" - -import collections -import ufl - -from gem import (Literal, Zero, Identity, Sum, Product, Division, - Power, MathFunction, MinValue, MaxValue, Comparison, - LogicalNot, LogicalAnd, LogicalOr, Conditional, - Index, Indexed, ComponentTensor, IndexSum, - ListTensor) - - -class Mixin(object): - """A mixin to be used with a UFL MultiFunction to translate UFL - algebra into GEM tensor-algebra. This node types translate pretty - straightforwardly to GEM. Other node types are not handled in - this mixin.""" - - def __init__(self): - self.index_map = collections.defaultdict(Index) - """A map for translating UFL free indices into GEM free - indices.""" - - def scalar_value(self, o): - return Literal(o.value()) - - def identity(self, o): - return Identity(o._dim) - - def zero(self, o): - return Zero(o.ufl_shape) - - def sum(self, o, *ops): - if o.ufl_shape: - indices = tuple(Index() for i in range(len(o.ufl_shape))) - return ComponentTensor(Sum(*[Indexed(op, indices) for op in ops]), indices) - else: - return Sum(*ops) - - def product(self, o, *ops): - assert o.ufl_shape == () - return Product(*ops) - - def division(self, o, numerator, denominator): - return Division(numerator, denominator) - - def real(self, o, expr): - if o.ufl_shape: - indices = tuple(Index() for i in range(len(o.ufl_shape))) - return ComponentTensor(MathFunction('real', Indexed(expr, indices)), indices) - else: - return MathFunction('real', expr) - - def imag(self, o, expr): - if o.ufl_shape: - indices = tuple(Index() for i in range(len(o.ufl_shape))) - return ComponentTensor(MathFunction('imag', Indexed(expr, indices)), indices) - else: - return MathFunction('imag', expr) - - def conj(self, o, expr): - if o.ufl_shape: - indices = tuple(Index() for i in range(len(o.ufl_shape))) - return ComponentTensor(MathFunction('conj', Indexed(expr, indices)), indices) - else: - return MathFunction('conj', expr) - - def abs(self, o, expr): - if o.ufl_shape: - indices = tuple(Index() for i in range(len(o.ufl_shape))) - return ComponentTensor(MathFunction('abs', Indexed(expr, indices)), indices) - else: - return MathFunction('abs', expr) - - def power(self, o, base, exponent): - return Power(base, exponent) - - def math_function(self, o, expr): - return MathFunction(o._name, expr) - - def atan2(self, o, y, x): - return MathFunction("atan2", y, x) - - def bessel_i(self, o, nu, arg): - return MathFunction(o._name, nu, arg) - - def bessel_j(self, o, nu, arg): - return MathFunction(o._name, nu, arg) - - def bessel_k(self, o, nu, arg): - return MathFunction(o._name, nu, arg) - - def bessel_y(self, o, nu, arg): - return MathFunction(o._name, nu, arg) - - def min_value(self, o, *ops): - return MinValue(*ops) - - def max_value(self, o, *ops): - return MaxValue(*ops) - - def binary_condition(self, o, left, right): - return Comparison(o._name, left, right) - - def not_condition(self, o, expr): - return LogicalNot(expr) - - def and_condition(self, o, *ops): - return LogicalAnd(*ops) - - def or_condition(self, o, *ops): - return LogicalOr(*ops) - - def conditional(self, o, condition, then, else_): - assert condition.shape == () - if o.ufl_shape: - indices = tuple(Index() for i in range(len(o.ufl_shape))) - return ComponentTensor(Conditional(condition, Indexed(then, indices), - Indexed(else_, indices)), - indices) - else: - return Conditional(condition, then, else_) - - def multi_index(self, o): - indices = [] - for i in o: - if isinstance(i, ufl.classes.FixedIndex): - indices.append(int(i)) - elif isinstance(i, ufl.classes.Index): - indices.append(self.index_map[i.count()]) - return tuple(indices) - - def indexed(self, o, aggregate, index): - return Indexed(aggregate, index) - - def list_tensor(self, o, *ops): - return ListTensor(ops) - - def component_tensor(self, o, expression, index): - return ComponentTensor(expression, index) - - def index_sum(self, o, summand, indices): - # ufl.IndexSum technically has a MultiIndex, but it must have - # exactly one index in it. - index, = indices - - if o.ufl_shape: - indices = tuple(Index() for i in range(len(o.ufl_shape))) - return ComponentTensor(IndexSum(Indexed(summand, indices), (index,)), indices) - else: - return IndexSum(summand, (index,)) - - def variable(self, o, expression, label): - # Only used by UFL AD, at this point, the bare expression is - # what we want. - return expression - - def label(self, o): - # Only used by UFL AD, don't need it at this point. - pass diff --git a/tsfc/ufl_utils.py b/tsfc/ufl_utils.py deleted file mode 100644 index 0cec4dbf..00000000 --- a/tsfc/ufl_utils.py +++ /dev/null @@ -1,485 +0,0 @@ -"""Utilities for preprocessing UFL objects.""" - -from functools import singledispatch - -import numpy - -import ufl -from ufl import as_tensor, indices, replace -from ufl.algorithms import compute_form_data as ufl_compute_form_data -from ufl.algorithms import estimate_total_polynomial_degree -from ufl.algorithms.analysis import extract_arguments, extract_type -from ufl.algorithms.apply_function_pullbacks import apply_function_pullbacks -from ufl.algorithms.apply_algebra_lowering import apply_algebra_lowering -from ufl.algorithms.apply_derivatives import apply_derivatives -from ufl.algorithms.apply_geometry_lowering import apply_geometry_lowering -from ufl.algorithms.apply_restrictions import apply_restrictions -from ufl.algorithms.comparison_checker import do_comparison_check -from ufl.algorithms.remove_complex_nodes import remove_complex_nodes -from ufl.corealg.map_dag import map_expr_dag -from ufl.corealg.multifunction import MultiFunction -from ufl.geometry import QuadratureWeight -from ufl.geometry import Jacobian, JacobianDeterminant, JacobianInverse -from ufl.classes import (Abs, Argument, CellOrientation, Coefficient, - ComponentTensor, Expr, FloatValue, Division, - MultiIndex, Product, - ScalarValue, Sqrt, Zero, CellVolume, FacetArea) -from ufl.domain import extract_unique_domain -from finat.ufl import MixedElement -from ufl.utils.sorting import sorted_by_count - -from gem.node import MemoizerArg - -from tsfc.modified_terminals import (is_modified_terminal, - analyse_modified_terminal, - construct_modified_terminal) - - -preserve_geometry_types = (CellVolume, FacetArea) - - -def compute_form_data(form, - do_apply_function_pullbacks=True, - do_apply_integral_scaling=True, - do_apply_geometry_lowering=True, - preserve_geometry_types=preserve_geometry_types, - do_apply_restrictions=True, - do_estimate_degrees=True, - complex_mode=False): - """Preprocess UFL form in a format suitable for TSFC. Return - form data. - - This is merely a wrapper to UFL compute_form_data with default - kwargs overriden in the way TSFC needs it and is provided for - other form compilers based on TSFC. - """ - fd = ufl_compute_form_data( - form, - do_apply_function_pullbacks=do_apply_function_pullbacks, - do_apply_integral_scaling=do_apply_integral_scaling, - do_apply_geometry_lowering=do_apply_geometry_lowering, - preserve_geometry_types=preserve_geometry_types, - do_apply_restrictions=do_apply_restrictions, - do_estimate_degrees=do_estimate_degrees, - complex_mode=complex_mode - ) - constants = extract_firedrake_constants(form) - fd.constants = constants - return fd - - -def extract_firedrake_constants(a): - """Build a sorted list of all constants in a""" - return sorted_by_count(extract_type(a, TSFCConstantMixin)) - - -def one_times(measure): - # Workaround for UFL issue #80: - # https://bitbucket.org/fenics-project/ufl/issues/80 - form = 1 * measure - fd = compute_form_data(form, do_estimate_degrees=False) - itg_data, = fd.integral_data - integral, = itg_data.integrals - integrand = integral.integrand() - - # UFL considers QuadratureWeight a geometric quantity, and the - # general handler for geometric quantities estimates the degree of - # the coordinate element. This would unnecessarily increase the - # estimated degree, so we drop QuadratureWeight instead. - expression = replace(integrand, {QuadratureWeight(itg_data.domain): 1}) - - # Now estimate degree for the preprocessed form - degree = estimate_total_polynomial_degree(expression) - - return integrand, degree - - -def entity_avg(integrand, measure, argument_multiindices): - arguments = extract_arguments(integrand) - if len(arguments) == 1: - a, = arguments - integrand = ufl.replace(integrand, {a: ufl.Argument(a.function_space(), - number=0, - part=a.part())}) - argument_multiindices = (argument_multiindices[a.number()], ) - - degree = estimate_total_polynomial_degree(integrand) - form = integrand * measure - fd = compute_form_data(form, do_estimate_degrees=False, - do_apply_function_pullbacks=False) - itg_data, = fd.integral_data - integral, = itg_data.integrals - integrand = integral.integrand() - return integrand, degree, argument_multiindices - - -def preprocess_expression(expression, complex_mode=False, - do_apply_restrictions=False): - """Imitates the compute_form_data processing pipeline. - - :arg complex_mode: Are we in complex UFL mode? - :arg do_apply_restrictions: Propogate restrictions to terminals? - - Useful, for example, to preprocess non-scalar expressions, which - are not and cannot be forms. - """ - if complex_mode: - expression = do_comparison_check(expression) - else: - expression = remove_complex_nodes(expression) - expression = apply_algebra_lowering(expression) - expression = apply_derivatives(expression) - expression = apply_function_pullbacks(expression) - expression = apply_geometry_lowering(expression, preserve_geometry_types) - expression = apply_derivatives(expression) - expression = apply_geometry_lowering(expression, preserve_geometry_types) - expression = apply_derivatives(expression) - if not complex_mode: - expression = remove_complex_nodes(expression) - if do_apply_restrictions: - expression = apply_restrictions(expression) - return expression - - -class ModifiedTerminalMixin(object): - """Mixin to use with MultiFunctions that operate on modified - terminals.""" - - def unexpected(self, o): - assert False, "Not expected %r at this stage." % o - - # global derivates should have been pulled back - grad = unexpected - div = unexpected - curl = unexpected - - # div and curl should have been algebraically lowered - reference_div = unexpected - reference_curl = unexpected - - def _modified_terminal(self, o): - assert is_modified_terminal(o) - return self.modified_terminal(o) - - # Unlike UFL, we do not regard Indexed as a terminal modifier. - # indexed = _modified_terminal - - positive_restricted = _modified_terminal - negative_restricted = _modified_terminal - - reference_grad = _modified_terminal - reference_value = _modified_terminal - - terminal = _modified_terminal - - -class CoefficientSplitter(MultiFunction, ModifiedTerminalMixin): - def __init__(self, split): - MultiFunction.__init__(self) - self._split = split - - expr = MultiFunction.reuse_if_untouched - - def modified_terminal(self, o): - mt = analyse_modified_terminal(o) - terminal = mt.terminal - - if not isinstance(terminal, Coefficient): - # Only split coefficients - return o - - if type(terminal.ufl_element()) != MixedElement: - # Only split mixed coefficients - return o - - # Reference value expected - assert mt.reference_value - - # Derivative indices - beta = indices(mt.local_derivatives) - - components = [] - for subcoeff in self._split[terminal]: - # Apply terminal modifiers onto the subcoefficient - component = construct_modified_terminal(mt, subcoeff) - # Collect components of the subcoefficient - for alpha in numpy.ndindex(subcoeff.ufl_element().reference_value_shape): - # New modified terminal: component[alpha + beta] - components.append(component[alpha + beta]) - # Repack derivative indices to shape - c, = indices(1) - return ComponentTensor(as_tensor(components)[c], MultiIndex((c,) + beta)) - - -def split_coefficients(expression, split): - """Split mixed coefficients, so mixed elements need not be - implemented. - - :arg split: A :py:class:`dict` mapping each mixed coefficient to a - sequence of subcoefficients. If None, calling this - function is a no-op. - """ - if split is None: - return expression - - splitter = CoefficientSplitter(split) - return map_expr_dag(splitter, expression) - - -class PickRestriction(MultiFunction, ModifiedTerminalMixin): - """Pick out parts of an expression with specified restrictions on - the arguments. - - :arg test: The restriction on the test function. - :arg trial: The restriction on the trial function. - - Returns those parts of the expression that have the requested - restrictions, or else :class:`ufl.classes.Zero` if no such part - exists. - """ - def __init__(self, test=None, trial=None): - self.restrictions = {0: test, 1: trial} - MultiFunction.__init__(self) - - expr = MultiFunction.reuse_if_untouched - - def multi_index(self, o): - return o - - def modified_terminal(self, o): - mt = analyse_modified_terminal(o) - t = mt.terminal - r = mt.restriction - if isinstance(t, Argument) and r != self.restrictions[t.number()]: - return Zero(o.ufl_shape, o.ufl_free_indices, o.ufl_index_dimensions) - else: - return o - - -def ufl_reuse_if_untouched(o, *ops): - """Reuse object if operands are the same objects.""" - if all(a is b for a, b in zip(o.ufl_operands, ops)): - return o - else: - return o._ufl_expr_reconstruct_(*ops) - - -@singledispatch -def _simplify_abs(o, self, in_abs): - """Single-dispatch function to simplify absolute values. - - :arg o: UFL node - :arg self: Callback handler for recursion - :arg in_abs: Is ``o`` inside an absolute value? - - When ``in_abs`` we must return a non-negative value, potentially - by wrapping the returned node with ``Abs``. - """ - raise AssertionError("UFL node expected, not %s" % type(o)) - - -@_simplify_abs.register(Expr) -def _simplify_abs_expr(o, self, in_abs): - # General case, only wrap the outer expression (if necessary) - operands = [self(op, False) for op in o.ufl_operands] - result = ufl_reuse_if_untouched(o, *operands) - if in_abs: - result = Abs(result) - return result - - -@_simplify_abs.register(Sqrt) -def _simplify_abs_sqrt(o, self, in_abs): - result = ufl_reuse_if_untouched(o, self(o.ufl_operands[0], False)) - if self.complex_mode and in_abs: - return Abs(result) - else: - return result - - -@_simplify_abs.register(ScalarValue) -def _simplify_abs_(o, self, in_abs): - if not in_abs: - return o - # Inline abs(constant) - return ufl.as_ufl(abs(o._value)) - - -@_simplify_abs.register(CellOrientation) -def _simplify_abs_cellorientation(o, self, in_abs): - if not in_abs: - return o - # Cell orientation is +-1 - return FloatValue(1) - - -@_simplify_abs.register(Division) -@_simplify_abs.register(Product) -def _simplify_abs_product(o, self, in_abs): - if not in_abs: - # Just reconstruct - ops = [self(op, False) for op in o.ufl_operands] - return ufl_reuse_if_untouched(o, *ops) - - # Visit children, distributing Abs - ops = [self(op, True) for op in o.ufl_operands] - - # Strip Abs off again (we will put it outside now) - stripped = False - strip_ops = [] - for op in ops: - if isinstance(op, Abs): - stripped = True - strip_ops.append(op.ufl_operands[0]) - else: - strip_ops.append(op) - - # Rebuild, and wrap with Abs if necessary - result = ufl_reuse_if_untouched(o, *strip_ops) - if stripped: - result = Abs(result) - return result - - -@_simplify_abs.register(Abs) -def _simplify_abs_abs(o, self, in_abs): - return self(o.ufl_operands[0], True) - - -def simplify_abs(expression, complex_mode): - """Simplify absolute values in a UFL expression. Its primary - purpose is to "neutralise" CellOrientation nodes that are - surrounded by absolute values and thus not at all necessary.""" - mapper = MemoizerArg(_simplify_abs) - mapper.complex_mode = complex_mode - return mapper(expression, False) - - -def apply_mapping(expression, element, domain): - """Apply the inverse of the pullback for element to an expression. - - :arg expression: An expression in physical space - :arg element: The element we're going to interpolate into, whose - value_shape must match the shape of the expression, and will - advertise the pullback to apply. - :arg domain: Optional domain to provide in case expression does - not contain a domain (used for constructing geometric quantities). - :returns: A new UFL expression with shape element.reference_value_shape - :raises NotImplementedError: If we don't know how to apply the - inverse of the pullback. - :raises ValueError: If we get shape mismatches. - - The following is borrowed from the UFC documentation: - - Let g be a field defined on a physical domain T with physical - coordinates x. Let T_0 be a reference domain with coordinates - X. Assume that F: T_0 -> T such that - - x = F(X) - - Let J be the Jacobian of F, i.e J = dx/dX and let K denote the - inverse of the Jacobian K = J^{-1}. Then we (currently) have the - following four types of mappings: - - 'identity' mapping for g: - - G(X) = g(x) - - For vector fields g: - - 'contravariant piola' mapping for g: - - G(X) = det(J) K g(x) i.e G_i(X) = det(J) K_ij g_j(x) - - 'covariant piola' mapping for g: - - G(X) = J^T g(x) i.e G_i(X) = J^T_ij g(x) = J_ji g_j(x) - - 'double covariant piola' mapping for g: - - G(X) = J^T g(x) J i.e. G_il(X) = J_ji g_jk(x) J_kl - - 'double contravariant piola' mapping for g: - - G(X) = det(J)^2 K g(x) K^T i.e. G_il(X)=(detJ)^2 K_ij g_jk K_lk - - If 'contravariant piola' or 'covariant piola' (or their double - variants) are applied to a matrix-valued function, the appropriate - mappings are applied row-by-row. - """ - mesh = extract_unique_domain(expression) - if mesh is None: - mesh = domain - if domain is not None and mesh != domain: - raise NotImplementedError("Multiple domains not supported") - if expression.ufl_shape != element.value_shape: - raise ValueError(f"Mismatching shapes, got {expression.ufl_shape}, expected {element.value_shape}") - mapping = element.mapping().lower() - if mapping == "identity": - rexpression = expression - elif mapping == "covariant piola": - J = Jacobian(mesh) - *k, i, j = indices(len(expression.ufl_shape) + 1) - kj = (*k, j) - rexpression = as_tensor(J[j, i] * expression[kj], (*k, i)) - elif mapping == "l2 piola": - detJ = JacobianDeterminant(mesh) - rexpression = expression * detJ - elif mapping == "contravariant piola": - K = JacobianInverse(mesh) - detJ = JacobianDeterminant(mesh) - *k, i, j = indices(len(expression.ufl_shape) + 1) - kj = (*k, j) - rexpression = as_tensor(detJ * K[i, j] * expression[kj], (*k, i)) - elif mapping == "double covariant piola": - J = Jacobian(mesh) - *k, i, j, m, n = indices(len(expression.ufl_shape) + 2) - kmn = (*k, m, n) - rexpression = as_tensor(J[m, i] * expression[kmn] * J[n, j], (*k, i, j)) - elif mapping == "double contravariant piola": - K = JacobianInverse(mesh) - detJ = JacobianDeterminant(mesh) - *k, i, j, m, n = indices(len(expression.ufl_shape) + 2) - kmn = (*k, m, n) - rexpression = as_tensor(detJ**2 * K[i, m] * expression[kmn] * K[j, n], (*k, i, j)) - elif mapping == "symmetries": - # This tells us how to get from the pieces of the reference - # space expression to the physical space one. - # We're going to apply the inverse of the physical to - # reference space mapping. - fcm = element.flattened_sub_element_mapping() - sub_elem = element.sub_elements[0] - shape = expression.ufl_shape - flat = ufl.as_vector([expression[i] for i in numpy.ndindex(shape)]) - vs = sub_elem.value_shape - rvs = sub_elem.reference_value_shape - seen = set() - rpieces = [] - gm = int(numpy.prod(vs, dtype=int)) - for gi, ri in enumerate(fcm): - # For each unique piece in reference space - if ri in seen: - continue - seen.add(ri) - # Get the physical space piece - piece = [flat[gm*gi + j] for j in range(gm)] - piece = as_tensor(numpy.asarray(piece).reshape(vs)) - # get into reference space - piece = apply_mapping(piece, sub_elem, mesh) - assert piece.ufl_shape == rvs - # Concatenate with the other pieces - rpieces.extend([piece[idx] for idx in numpy.ndindex(rvs)]) - # And reshape - rexpression = as_tensor(numpy.asarray(rpieces).reshape(element.reference_value_shape)) - else: - raise NotImplementedError(f"Don't know how to handle mapping type {mapping} for expression of rank {element.value_shape}") - if rexpression.ufl_shape != element.reference_value_shape: - raise ValueError(f"Mismatching reference shapes, got {rexpression.ufl_shape} expected {element.reference_value_shape}") - return rexpression - - -class TSFCConstantMixin: - """ Mixin class to identify Constants """ - - def __init__(self): - pass diff --git a/tsfc/vanilla.py b/tsfc/vanilla.py deleted file mode 100644 index ecf24cd5..00000000 --- a/tsfc/vanilla.py +++ /dev/null @@ -1,47 +0,0 @@ -from functools import reduce - -from gem import index_sum, Sum -from gem.optimise import unroll_indexsum -from gem.unconcatenate import unconcatenate - - -def Integrals(expressions, quadrature_multiindex, argument_multiindices, parameters): - """Constructs an integral representation for each GEM integrand - expression. - - :arg expressions: integrand multiplied with quadrature weight; - multi-root GEM expression DAG - :arg quadrature_multiindex: quadrature multiindex (tuple) - :arg argument_multiindices: tuple of argument multiindices, - one multiindex for each argument - :arg parameters: parameters dictionary - - :returns: list of integral representations - """ - # Unroll - max_extent = parameters["unroll_indexsum"] - if max_extent: - def predicate(index): - return index.extent <= max_extent - expressions = unroll_indexsum(expressions, predicate=predicate) - # Integral representation: just a GEM expression - return [index_sum(e, quadrature_multiindex) for e in expressions] - - -def flatten(var_reps, index_cache): - """Flatten mode-specific intermediate representation to a series of - assignments. - - :arg var_reps: series of (return variable, [integral representation]) pairs - :arg index_cache: cache :py:class:`dict` for :py:func:`unconcatenate` - - :returns: series of (return variable, GEM expression root) pairs - """ - return unconcatenate([(variable, reduce(Sum, reps)) - for variable, reps in var_reps], - cache=index_cache) - - -finalise_options = dict(remove_componenttensors=False) -"""To avoid duplicate work, these options that are safe to pass to -:py:func:`gem.impero_utils.preprocess_gem`."""