From b39f7f6082a11444abe482e191f7e12d2bae309a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 09:51:24 +0000 Subject: [PATCH 01/10] Introduce integral type and try to separate ufl concepts from FFCx concepts --- ffcx/analysis.py | 6 +- ffcx/codegeneration/C/form.py | 11 +- ffcx/codegeneration/access.py | 70 +++++------ ffcx/codegeneration/backend.py | 8 +- ffcx/codegeneration/definitions.py | 17 ++- ffcx/codegeneration/expression_generator.py | 4 +- ffcx/codegeneration/integral_generator.py | 3 +- ffcx/codegeneration/jit.py | 2 +- ffcx/codegeneration/symbols.py | 44 +++---- ffcx/codegeneration/ufcx.h | 2 +- ffcx/definitions.py | 38 +++++- ffcx/ir/elementtables.py | 101 ++++++++------- ffcx/ir/integral.py | 11 +- ffcx/ir/representation.py | 88 ++++++------- ffcx/ir/representationutils.py | 130 ++++++++++---------- test/test_add_mode.py | 2 +- test/test_jit_forms.py | 6 +- 17 files changed, 290 insertions(+), 253 deletions(-) diff --git a/ffcx/analysis.py b/ffcx/analysis.py index 31d214012..dba9c6be3 100644 --- a/ffcx/analysis.py +++ b/ffcx/analysis.py @@ -93,8 +93,8 @@ def analyze_ufl_objects( form_data = tuple(_analyze_form(form, scalar_type) for form in forms) for data in form_data: - elements += data.unique_sub_elements # type: ignore - coordinate_elements += data.coordinate_elements # type: ignore + elements += data.unique_sub_elements + coordinate_elements += data.coordinate_elements for original_expression, points in expressions: elements += ufl.algorithms.extract_elements(original_expression) @@ -190,7 +190,7 @@ def _analyze_form( # Determine unique quadrature degree and quadrature scheme # per each integral data - for id, integral_data in enumerate(form_data.integral_data): # type: ignore + for id, integral_data in enumerate(form_data.integral_data): # Iterate through groups of integral data. There is one integral # data for all integrals with same domain, itype, subdomain_id # (but possibly different metadata). diff --git a/ffcx/codegeneration/C/form.py b/ffcx/codegeneration/C/form.py index f026dc386..8d33a4044 100644 --- a/ffcx/codegeneration/C/form.py +++ b/ffcx/codegeneration/C/form.py @@ -4,7 +4,7 @@ # # SPDX-License-Identifier: LGPL-3.0-or-later # -# Modified by Chris Richardson and Jørgen S. Dokken 2023 +# Modified by Chris Richardson and Jørgen S. Dokken 2023-2025 # # Note: Most of the code in this file is a direct translation from the # old implementation in FFC @@ -17,6 +17,7 @@ import numpy as np from ffcx.codegeneration.C import form_template +from ffcx.definitions import IntegralType from ffcx.ir.representation import FormIR logger = logging.getLogger("ffcx") @@ -119,7 +120,13 @@ def generator(ir: FormIR, options): integral_offsets = [0] integral_domains = [] # Note: the order of this list is defined by the enum ufcx_integral_type in ufcx.h - for itg_type in ("cell", "exterior_facet", "interior_facet", "vertex", "ridge"): + for itg_type in ( + IntegralType(codim=0), + IntegralType(codim=1, num_neighbours=1), + IntegralType(codim=1, num_neighbours=2), + IntegralType(codim=-1), + IntegralType(codim=2), + ): unsorted_integrals = [] unsorted_ids = [] unsorted_domains = [] diff --git a/ffcx/codegeneration/access.py b/ffcx/codegeneration/access.py index 78ddb5566..3125a71ef 100644 --- a/ffcx/codegeneration/access.py +++ b/ffcx/codegeneration/access.py @@ -1,4 +1,4 @@ -# Copyright (C) 2011-2017 Martin Sandve Alnæs +# Copyright (C) 2011-2025 Martin Sandve Alnæs, Jørgen S. Dokken # # This file is part of FFCx. (https://www.fenicsproject.org) # @@ -12,7 +12,7 @@ import ufl import ffcx.codegeneration.lnodes as L -from ffcx.definitions import entity_types +from ffcx.definitions import IntegralType from ffcx.ir.analysis.modified_terminals import ModifiedTerminal from ffcx.ir.elementtables import UniqueTableReferenceT from ffcx.ir.representationutils import QuadratureRule @@ -23,12 +23,11 @@ class FFCXBackendAccess: """FFCx specific formatter class.""" - entity_type: entity_types + integral_type: IntegralType - def __init__(self, entity_type: entity_types, integral_type: str, symbols, options): + def __init__(self, integral_type: IntegralType, symbols, options): """Initialise.""" # Store ir and options - self.entity_type = entity_type self.integral_type = integral_type self.symbols = symbols self.options = options @@ -65,7 +64,7 @@ def get( """Format a terminal.""" e = mt.terminal # Call appropriate handler, depending on the type of e - handler = self.call_lookup.get(type(e), False) # type: ignore + handler = self.call_lookup.get(type(e), False) if not handler: # Look for parent class types instead @@ -74,7 +73,7 @@ def get( handler = self.call_lookup[k] break if handler: - return handler(mt, tabledata, quadrature_rule) + return handler(mt, tabledata, quadrature_rule) # type: ignore else: raise RuntimeError(f"Not handled: {type(e)}") @@ -122,27 +121,28 @@ def spatial_coordinate( if mt.averaged is not None: raise RuntimeError("Not expecting average of SpatialCoordinates.") - if self.integral_type in ufl.custom_integral_types: - if mt.local_derivatives: - raise RuntimeError("FIXME: Jacobian in custom integrals is not implemented.") - - # Access predefined quadrature points table - x = self.symbols.custom_points_table - iq = self.symbols.quadrature_loop_index - (gdim,) = mt.terminal.ufl_shape - if gdim == 1: - index = iq - else: - index = iq * gdim + mt.flat_component - return x[index] - elif self.integral_type == "expression": - # Physical coordinates are computed by code generated in - # definitions - return self.symbols.x_component(mt) - else: - # Physical coordinates are computed by code generated in - # definitions - return self.symbols.x_component(mt) + match self.integral_type: + case IntegralType(is_custom=True): + if mt.local_derivatives: + raise RuntimeError("FIXME: Jacobian in custom integrals is not implemented.") + + # Access predefined quadrature points table + x = self.symbols.custom_points_table + iq = self.symbols.quadrature_loop_index + (gdim,) = mt.terminal.ufl_shape + if gdim == 1: + index = iq + else: + index = iq * gdim + mt.flat_component + return x[index] + case IntegralType(is_expression=True): + # Physical coordinates are computed by code generated in + # definitions + return self.symbols.x_component(mt) + case _: + # Physical coordinates are computed by code generated in + # definitions + return self.symbols.x_component(mt) def cell_coordinate(self, mt, tabledata, num_points): """Access a cell coordinate.""" @@ -181,7 +181,7 @@ def facet_coordinate(self, mt, tabledata, num_points): if mt.restriction: raise RuntimeError("Not expecting restriction of FacetCoordinate.") - if self.integral_type in ("interior_facet", "exterior_facet"): + if self.integral_type in ("interior_facet", "facet"): (tdim,) = mt.terminal.ufl_shape if tdim == 0: raise RuntimeError("Vertices have no facet coordinates.") @@ -233,7 +233,7 @@ def reference_normal(self, mt, tabledata, access): cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname() if cellname in ("interval", "triangle", "tetrahedron", "quadrilateral", "hexahedron"): table = L.Symbol(f"{cellname}_reference_normals", dtype=L.DataType.REAL) - facet = self.symbols.entity("facet", mt.restriction) + facet = self.symbols.entity(IntegralType(codim=1), mt.restriction) return table[facet][mt.component[0]] else: raise RuntimeError(f"Unhandled cell types {cellname}.") @@ -250,7 +250,7 @@ def cell_facet_jacobian(self, mt, tabledata, num_points): "pyramid", ): table = L.Symbol(f"{cellname}_cell_facet_jacobian", dtype=L.DataType.REAL) - facet = self.symbols.entity("facet", mt.restriction) + facet = self.symbols.entity(IntegralType(codim=1), mt.restriction) return table[facet][mt.component[0]][mt.component[1]] elif cellname == "interval": raise RuntimeError("The reference facet jacobian doesn't make sense for interval cell.") @@ -262,7 +262,7 @@ def cell_ridge_jacobian(self, mt, tabledata, num_points): cellname = ufl.domain.extract_unique_domain(mt.terminal).ufl_cell().cellname() if cellname in ("tetrahedron", "prism", "hexahedron"): table = L.Symbol(f"{cellname}_cell_ridge_jacobian", dtype=L.DataType.REAL) - ridge = self.symbols.entity("ridge", mt.restriction) + ridge = self.symbols.entity(IntegralType(codim=2), mt.restriction) return table[ridge][mt.component[0]][mt.component[1]] elif cellname in ["triangle", "quadrilateral"]: raise RuntimeError("The ridge jacobian doesn't make sense for 2D cells.") @@ -422,7 +422,7 @@ def _pass(self, *args, **kwargs): def table_access( self, tabledata: UniqueTableReferenceT, - entity_type: entity_types, + integral_type: IntegralType, restriction: str, quadrature_index: L.MultiIndex, dof_index: L.MultiIndex, @@ -431,12 +431,12 @@ def table_access( Args: tabledata: Table data object - entity_type: Entity type + integral_type: Integral type restriction: Restriction ("+", "-") quadrature_index: Quadrature index dof_index: Dof index """ - entity = self.symbols.entity(entity_type, restriction) + entity = self.symbols.entity(integral_type, restriction) iq_global_index = quadrature_index.global_index ic_global_index = dof_index.global_index qp = 0 # quadrature permutation diff --git a/ffcx/codegeneration/backend.py b/ffcx/codegeneration/backend.py index ef6963987..165dcf4cb 100644 --- a/ffcx/codegeneration/backend.py +++ b/ffcx/codegeneration/backend.py @@ -26,9 +26,5 @@ def __init__(self, ir: IntegralIR | ExpressionIR, options): self.symbols = FFCXBackendSymbols( coefficient_numbering, coefficient_offsets, original_constant_offsets ) - self.access = FFCXBackendAccess( - ir.expression.entity_type, ir.expression.integral_type, self.symbols, options - ) - self.definitions = FFCXBackendDefinitions( - ir.expression.entity_type, ir.expression.integral_type, self.access, options - ) + self.access = FFCXBackendAccess(ir.expression.integral_type, self.symbols, options) + self.definitions = FFCXBackendDefinitions(ir.expression.integral_type, self.access, options) diff --git a/ffcx/codegeneration/definitions.py b/ffcx/codegeneration/definitions.py index dd207db4c..3f4b13220 100644 --- a/ffcx/codegeneration/definitions.py +++ b/ffcx/codegeneration/definitions.py @@ -1,4 +1,4 @@ -# Copyright (C) 2011-2023 Martin Sandve Alnæs, Igor A. Baratta +# Copyright (C) 2011-2025 Martin Sandve Alnæs, Igor A. Baratta and Jørgen S. Dokken # # This file is part of FFCx. (https://www.fenicsproject.org) # @@ -10,7 +10,7 @@ import ufl import ffcx.codegeneration.lnodes as L -from ffcx.definitions import entity_types +from ffcx.definitions import IntegralType from ffcx.ir.analysis.modified_terminals import ModifiedTerminal from ffcx.ir.elementtables import UniqueTableReferenceT from ffcx.ir.representationutils import QuadratureRule @@ -50,13 +50,12 @@ def create_dof_index(tabledata, dof_index_symbol): class FFCXBackendDefinitions: """FFCx specific code definitions.""" - entity_type: entity_types + integral_type: IntegralType - def __init__(self, entity_type: entity_types, integral_type: str, access, options): + def __init__(self, integral_type: IntegralType, access, options): """Initialise.""" # Store ir and options self.integral_type = integral_type - self.entity_type = entity_type self.access = access self.options = options @@ -102,13 +101,13 @@ def get( ttype = ttype.__bases__[0] # Get the handler from the lookup, or None if not found - handler = self.handler_lookup.get(ttype) # type: ignore + handler = self.handler_lookup.get(ttype) if handler is None: raise NotImplementedError(f"No handler for terminal type: {ttype}") # Call the handler - return handler(mt, tabledata, quadrature_rule, access) # type: ignore + return handler(mt, tabledata, quadrature_rule, access) def coefficient( self, @@ -148,7 +147,7 @@ def coefficient( assert begin < end # Get access to element table - FE, tables = self.access.table_access(tabledata, self.entity_type, mt.restriction, iq, ic) + FE, tables = self.access.table_access(tabledata, self.integral_type, mt.restriction, iq, ic) dof_access: L.ArrayAccess = self.symbols.coefficient_dof_access( mt.terminal, (ic.global_index) * bs + begin ) @@ -197,7 +196,7 @@ def _define_coordinate_dofs_lincomb( iq_symbol = self.symbols.quadrature_loop_index ic = create_dof_index(tabledata, ic_symbol) iq = create_quadrature_index(quadrature_rule, iq_symbol) - FE, tables = self.access.table_access(tabledata, self.entity_type, mt.restriction, iq, ic) + FE, tables = self.access.table_access(tabledata, self.integral_type, mt.restriction, iq, ic) dof_access = L.Symbol("coordinate_dofs", dtype=L.DataType.REAL) diff --git a/ffcx/codegeneration/expression_generator.py b/ffcx/codegeneration/expression_generator.py index beab0fe45..bff079b3b 100644 --- a/ffcx/codegeneration/expression_generator.py +++ b/ffcx/codegeneration/expression_generator.py @@ -77,7 +77,7 @@ def generate_geometry_tables(self): mt = attr.get("mt") if mt is not None: t = type(mt.terminal) - if self.ir.expression.entity_type == "cell" and issubclass( + if self.ir.expression.integral_type.codim == 0 and issubclass( t, ufl.geometry.GeometricFacetQuantity ): raise RuntimeError(f"Expressions for cells do not support {t}.") @@ -302,7 +302,7 @@ def get_arg_factors(self, blockdata, block_rank, indices): ] table = self.backend.symbols.element_table( - td, self.ir.expression.entity_type, mt.restriction + td, self.ir.expression.integral_type, mt.restriction ) assert td.ttype != "zeros" diff --git a/ffcx/codegeneration/integral_generator.py b/ffcx/codegeneration/integral_generator.py index 45a3414cc..d7163c421 100644 --- a/ffcx/codegeneration/integral_generator.py +++ b/ffcx/codegeneration/integral_generator.py @@ -417,13 +417,12 @@ def get_arg_factors(self, blockdata, block_rank, quadrature_rule, domain, iq, in # now because it assumes too much about indices. assert td.ttype != "zeros" - if td.ttype == "ones": arg_factor = 1 else: # Assuming B sparsity follows element table sparsity arg_factor, arg_tables = self.backend.access.table_access( - td, self.ir.expression.entity_type, mt.restriction, iq, indices[i] + td, self.ir.expression.integral_type, mt.restriction, iq, indices[i] ) tables += arg_tables diff --git a/ffcx/codegeneration/jit.py b/ffcx/codegeneration/jit.py index cfedb5b3f..6eb5dbb8f 100644 --- a/ffcx/codegeneration/jit.py +++ b/ffcx/codegeneration/jit.py @@ -229,7 +229,7 @@ def compile_forms( def compile_expressions( - expressions: list[tuple[ufl.Expr, npt.NDArray[np.floating]]], # type: ignore + expressions: list[tuple[ufl.Expr, npt.NDArray[np.floating]]], options: dict = {}, cache_dir: Path | None = None, timeout: int = 10, diff --git a/ffcx/codegeneration/symbols.py b/ffcx/codegeneration/symbols.py index 4fb098a4c..9f0990f77 100644 --- a/ffcx/codegeneration/symbols.py +++ b/ffcx/codegeneration/symbols.py @@ -1,4 +1,4 @@ -# Copyright (C) 2011-2023 Martin Sandve Alnæs, Igor A. Baratta +# Copyright (C) 2011-2025 Martin Sandve Alnæs, Igor A. Baratta and Jørgen S. Dokken # # This file is part of FFCx. (https://www.fenicsproject.org) # @@ -6,11 +6,12 @@ """FFCx/UFC specific symbol naming.""" import logging +import typing import ufl import ffcx.codegeneration.lnodes as L -from ffcx.definitions import entity_types +from ffcx.definitions import IntegralType logger = logging.getLogger("ffcx") @@ -96,23 +97,25 @@ def __init__(self, coefficient_numbering, coefficient_offsets, original_constant # Table for chunk of custom quadrature points (physical coordinates). self.custom_points_table = L.Symbol("points_chunk", dtype=L.DataType.REAL) - def entity(self, entity_type: entity_types, restriction): + def entity(self, integral_type: IntegralType, restriction: typing.Literal["+", "-"]): """Entity index for lookup in element tables.""" - if entity_type == "cell": - # Always 0 for cells (even with restriction) - return L.LiteralInt(0) - - if entity_type == "facet": - if restriction == "-": - return self.entity_local_index[1] - else: + match integral_type: + case IntegralType(codim=0): + # Always 0 for cells (even with restriction) + return L.LiteralInt(0) + case IntegralType(codim=1): + if restriction == "-": + return self.entity_local_index[1] + else: + return self.entity_local_index[0] + case IntegralType(codim=-1): return self.entity_local_index[0] - elif entity_type == "vertex": - return self.entity_local_index[0] - elif entity_type == "ridge": - return self.entity_local_index[0] - else: - logger.exception(f"Unknown entity_type {entity_type}") + case IntegralType(codim=2): + return self.entity_local_index[0] + case _: + raise RuntimeError( + f"Unknown integral_type {integral_type} in entity restriction lookup/" + ) def argument_loop_index(self, iarg): """Loop index for argument iarg.""" @@ -180,14 +183,13 @@ def constant_index_access(self, constant, index): return c[offset + index] # TODO: Remove this, use table_access instead - def element_table(self, tabledata, entity_type: entity_types, restriction): + def element_table(self, tabledata, integral_type: IntegralType, restriction): """Get an element table.""" - entity = self.entity(entity_type, restriction) - + entity = self.entity(integral_type, restriction) if tabledata.is_uniform: entity = 0 else: - entity = self.entity(entity_type, restriction) + entity = self.entity(integral_type, restriction) if tabledata.is_piecewise: iq = 0 diff --git a/ffcx/codegeneration/ufcx.h b/ffcx/codegeneration/ufcx.h index 453cc9f89..9c85007a4 100644 --- a/ffcx/codegeneration/ufcx.h +++ b/ffcx/codegeneration/ufcx.h @@ -48,7 +48,7 @@ extern "C" typedef enum { cell = 0, - exterior_facet = 1, + facet = 1, interior_facet = 2, vertex = 3, ridge = 4, diff --git a/ffcx/definitions.py b/ffcx/definitions.py index 6aa30661f..7c304e971 100644 --- a/ffcx/definitions.py +++ b/ffcx/definitions.py @@ -1,5 +1,39 @@ +# Copyright (C) 2025 Jørgen S. Dokken +# +# This file is part of FFCx. (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later + """Module for storing type definitions used in the FFCx code base.""" -from typing import Literal +from typing import NamedTuple + + +class IntegralType(NamedTuple): + """Class for storing information about an integral type.""" + + codim: int # Codimension of integral with respect to the topological dimension of the domain + num_neighbours: int = 1 # Number of neighbouring cells + is_expression: bool = ( + False # True if the integral is an Expression (integrand) rather than an integral + ) + is_custom: bool = False # True if the integral uses the custom ufl integral measure. + -entity_types = Literal["cell", "facet", "vertex", "ridge"] +def convert_to_integral_type(integral_type: str) -> IntegralType: + """Convert UFL integral type string to IntegralType.""" + match integral_type: + case "cell": + return IntegralType(codim=0, num_neighbours=1) + case "exterior_facet": + return IntegralType(codim=1, num_neighbours=1) + case "interior_facet": + return IntegralType(codim=1, num_neighbours=2) + case "ridge": + return IntegralType(codim=2, num_neighbours=1) + case "vertex": + return IntegralType(codim=-1, num_neighbours=1) + case "custom": + return IntegralType(codim=0, num_neighbours=1, is_custom=True) + case _: + raise ValueError(f"Unknown integral type:{integral_type}") diff --git a/ffcx/ir/elementtables.py b/ffcx/ir/elementtables.py index 33698a506..e6582f8c9 100644 --- a/ffcx/ir/elementtables.py +++ b/ffcx/ir/elementtables.py @@ -1,4 +1,4 @@ -# Copyright (C) 2013-2017 Martin Sandve Alnæs +# Copyright (C) 2013-2025 Martin Sandve Alnæs and Jørgen S. Dokken # # This file is part of FFCx. (https://www.fenicsproject.org) # @@ -13,7 +13,7 @@ import numpy.typing as npt import ufl -from ffcx.definitions import entity_types +from ffcx.definitions import IntegralType from ffcx.element_interface import basix_index from ffcx.ir.analysis.modified_terminals import ModifiedTerminal from ffcx.ir.representationutils import ( @@ -80,12 +80,11 @@ def clamp_table_small_numbers( def get_ffcx_table_values( - points, - cell, - integral_type, - element, - avg, - entity_type: entity_types, + points: npt.NDArray[np.floating], + cell: ufl.Cell, + integral_type: IntegralType, + element: basix.ufl._ElementBase, + avg: typing.Literal["cell", "facet", None], derivative_counts, flat_component, codim, @@ -97,18 +96,14 @@ def get_ffcx_table_values( """ deriv_order = sum(derivative_counts) - if integral_type in ufl.custom_integral_types: - # Use quadrature points on cell for analysis in custom integral types - integral_type = "cell" - assert not avg - - if integral_type == "expression": - # FFCx tables for expression are generated as either interior cell points - # or points on a facet - if entity_type == "cell": - integral_type = "cell" - else: - integral_type = "exterior_facet" + match integral_type: + case IntegralType(is_custom=True): + assert integral_type.codim == 0 + integral_type = IntegralType(codim=0) + assert not avg + case IntegralType(is_expression=True): + assert integral_type.codim in (0, 1) + assert integral_type.num_neighbours == 1 if avg in ("cell", "facet"): # Redefine points to compute average tables @@ -123,17 +118,16 @@ def get_ffcx_table_values( # Doesn't matter if it's exterior or interior facet integral, # just need a valid integral type to create quadrature rule if avg == "cell": - integral_type = "cell" + integral_type = IntegralType(codim=0) elif avg == "facet": - integral_type = "exterior_facet" - + integral_type = IntegralType(codim=1) if isinstance(element, basix.ufl._QuadratureElement): points = element._points weights = element._weights else: # Make quadrature rule and get points and weights points, weights = create_quadrature_points_and_weights( - integral_type, cell, element.embedded_superdegree(), "default", [element] + integral_type, cell, element.embedded_superdegree, "default", [element] ) # Tabulate table of basis functions and derivatives in points for each entity @@ -179,8 +173,8 @@ def get_ffcx_table_values( def generate_psi_table_name( quadrature_rule: QuadratureRule, element_counter, - averaged: str, - entity_type: entity_types, + averaged: typing.Literal["cell", "facet", None], + integral_type: IntegralType, derivative_counts, flat_component, ): @@ -206,8 +200,21 @@ def generate_psi_table_name( name += f"_C{flat_component:d}" if any(derivative_counts): name += "_D" + "".join(str(d) for d in derivative_counts) - name += {None: "", "cell": "_AC", "facet": "_AF"}[averaged] - name += {"cell": "", "facet": "_F", "vertex": "_V", "ridge": "_R"}[entity_type] + + match averaged: + case "cell": + name += "_AC" + case "facet": + name += "_AF" + + match integral_type: + case IntegralType(codim=1): + name += "_F" + case IntegralType(codim=2): + name += "_R" + case IntegralType(codim=-1): + name += "_V" + name += f"_Q{quadrature_rule.id()}" return name @@ -223,7 +230,7 @@ def get_modified_terminal_element(mt) -> ModifiedTerminalElement | None: raise RuntimeError("Global derivatives of reference values not defined.") elif ld and not mt.reference_value: raise RuntimeError("Local derivatives of global values not defined.") - element = mt.terminal.ufl_function_space().ufl_element() # type: ignore + element = mt.terminal.ufl_function_space().ufl_element() fc = mt.flat_component elif isinstance(mt.terminal, ufl.classes.SpatialCoordinate): if mt.reference_value: @@ -315,8 +322,7 @@ def permute_quadrature_quadrilateral(points, reflections=0, rotations=0): def build_optimized_tables( quadrature_rule: QuadratureRule, cell: ufl.Cell, - integral_type: typing.Literal["interior_facet", "exterior_facet", "ridge", "cell", "vertex"], - entity_type: entity_types, + integral_type: IntegralType, modified_terminals: typing.Iterable[ModifiedTerminal], existing_tables: dict[str, npt.NDArray[np.float64]], use_sum_factorization: bool, @@ -329,7 +335,6 @@ def build_optimized_tables( Args: quadrature_rule: The quadrature rule used for the tables. cell: The cell type of the domain the tables will be used with. - entity_type: The entity type (vertex,edge,facet,cell) that the tables are evaluated for. integral_type: The type of integral the tables are used for. modified_terminals: Ordered sequence of unique modified terminals existing_tables: Register of tables that already exist and reused. @@ -375,8 +380,9 @@ def build_optimized_tables( # Build name for this particular table element_number = element_numbers[element] + name = generate_psi_table_name( - quadrature_rule, element_number, avg, entity_type, local_derivatives, flat_component + quadrature_rule, element_number, avg, integral_type, local_derivatives, flat_component # type: ignore ) # FIXME - currently just recalculate the tables every time, @@ -395,11 +401,11 @@ def build_optimized_tables( # needed because a cell may see its sub-entities as being oriented # differently to their global orientation if ( - integral_type == "interior_facet" - or integral_type == "ridge" + (integral_type.num_neighbours > 1 and integral_type.codim == 1) + or integral_type.codim == 2 or (is_mixed_dim and codim == 0) ): - if entity_type == "facet": + if integral_type.codim == 1: if tdim == 1 or codim == 1: # Do not add permutations if codim-1 as facets have already gotten a global # orientation in DOLFINx @@ -408,8 +414,7 @@ def build_optimized_tables( cell, integral_type, element, - avg, - entity_type, + avg, # type: ignore local_derivatives, flat_component, codim, @@ -423,8 +428,7 @@ def build_optimized_tables( cell, integral_type, element, - avg, - entity_type, + avg, # type: ignore local_derivatives, flat_component, codim, @@ -447,8 +451,7 @@ def build_optimized_tables( cell, integral_type, element, - avg, - entity_type, + avg, # type: ignore local_derivatives, flat_component, codim, @@ -468,8 +471,7 @@ def build_optimized_tables( cell, integral_type, element, - avg, - entity_type, + avg, # type: ignore local_derivatives, flat_component, codim, @@ -477,7 +479,7 @@ def build_optimized_tables( ) t = new_table[0] t["array"] = np.vstack([td["array"] for td in new_table]) - elif entity_type == "ridge": + elif integral_type.codim == 2: if tdim < 3 or codim == 2: # If ridge integral over vertex no permutation is needed, # or if it is a single domain ridge integral, @@ -487,8 +489,7 @@ def build_optimized_tables( cell, integral_type, element, - avg, - entity_type, + avg, # type: ignore local_derivatives, flat_component, codim, @@ -502,8 +503,7 @@ def build_optimized_tables( cell, integral_type, element, - avg, - entity_type, + avg, # type: ignore local_derivatives, flat_component, codim, @@ -517,8 +517,7 @@ def build_optimized_tables( cell, integral_type, element, - avg, - entity_type, + avg, # type: ignore local_derivatives, flat_component, codim, diff --git a/ffcx/ir/integral.py b/ffcx/ir/integral.py index 291badbc9..1fc800e59 100644 --- a/ffcx/ir/integral.py +++ b/ffcx/ir/integral.py @@ -18,7 +18,7 @@ from ufl.checks import is_cellwise_constant from ufl.classes import QuadratureWeight -from ffcx.definitions import entity_types +from ffcx.definitions import IntegralType from ffcx.ir.analysis.factorization import compute_argument_factorization from ffcx.ir.analysis.graph import build_scalar_graph from ffcx.ir.analysis.modified_terminals import ( @@ -36,8 +36,7 @@ class CommonExpressionIR(typing.NamedTuple): """Common-ground for IntegralIR and ExpressionIR.""" - integral_type: str - entity_type: entity_types + integral_type: IntegralType tensor_shape: list[int] coefficient_numbering: dict[ufl.Coefficient, int] coefficient_offsets: dict[ufl.Coefficient, int] @@ -74,8 +73,7 @@ class BlockDataT(typing.NamedTuple): def compute_integral_ir( cell: ufl.Cell, - integral_type: typing.Literal["interior_facet", "exterior_facet", "ridge", "cell"], - entity_type: typing.Literal["cell", "facet", "ridge", "vertex"], + integral_type: IntegralType, integrands: dict[basix.CellType, dict[QuadratureRule, ufl.core.expr.Expr]], argument_shape: tuple[int], p: dict, @@ -86,7 +84,6 @@ def compute_integral_ir( Args: cell: Cell of integration domain integral_type: Type of integral over cell - entity_type: Corresponding entity of the cell that the integral is over integrands: Dictionary mapping a quadrature rule to a sequence of integrands argument_shape: Shape of the output tensor of the integral (used for tensor factorization) p: Parameters used for clamping tables and for activating sum factorization @@ -135,12 +132,10 @@ def compute_integral_ir( for domain in ufl.domain.extract_domains(integrand): if domain.topological_dimension() != cell.topological_dimension(): is_mixed_dim = True - mt_table_reference = build_optimized_tables( quadrature_rule, cell, integral_type, - entity_type, initial_terminals.values(), ir["unique_tables"][integral_domain], use_sum_factorization=p["sum_factorization"], diff --git a/ffcx/ir/representation.py b/ffcx/ir/representation.py index fcd73c4d4..f5fe77d82 100644 --- a/ffcx/ir/representation.py +++ b/ffcx/ir/representation.py @@ -1,4 +1,4 @@ -# Copyright (C) 2009-2020 Anders Logg, Martin Sandve Alnæs, Marie E. Rognes, +# Copyright (C) 2009-2025 Anders Logg, Martin Sandve Alnæs, Marie E. Rognes, # Kristian B. Oelgaard, Matthew W. Scroggs, Chris Richardson, and others # # This file is part of FFCx. (https://www.fenicsproject.org) @@ -32,6 +32,7 @@ from ffcx import naming from ffcx.analysis import UFLData +from ffcx.definitions import IntegralType, convert_to_integral_type from ffcx.ir.integral import CommonExpressionIR, compute_integral_ir from ffcx.ir.representationutils import QuadratureRule, create_quadrature_points_and_weights @@ -78,9 +79,9 @@ class FormIR(typing.NamedTuple): constant_shapes: list[list[int]] constant_names: list[str] finite_element_hashes: list[int] - integral_names: dict[str, list[str]] - integral_domains: dict[str, list[basix.CellType]] - subdomain_ids: dict[str, list[int]] + integral_names: dict[IntegralType, list[str]] + integral_domains: dict[IntegralType, list[basix.CellType]] + subdomain_ids: dict[IntegralType, list[int]] class QuadratureIR(typing.NamedTuple): @@ -136,16 +137,15 @@ def compute_ir( integral_names = {} form_names = {} for fd_index, fd in enumerate(analysis.form_data): - form_names[fd_index] = naming.form_name(fd.original_form, fd_index, prefix) # type: ignore - for itg_index, itg_data in enumerate(fd.integral_data): # type: ignore + form_names[fd_index] = naming.form_name(fd.original_form, fd_index, prefix) + for itg_index, itg_data in enumerate(fd.integral_data): integral_names[(fd_index, itg_index)] = naming.integral_name( - fd.original_form, # type: ignore + fd.original_form, itg_data.integral_type, fd_index, itg_data.subdomain_id, prefix, ) - irs = [ _compute_integral_ir( fd, @@ -205,15 +205,6 @@ def _compute_integral_ir( visualise, ) -> list[IntegralIR]: """Compute intermediate representation for form integrals.""" - _entity_types = { - "cell": "cell", - "exterior_facet": "facet", - "interior_facet": "facet", - "vertex": "vertex", - "custom": "cell", - "ridge": "ridge", - } - # Iterate over groups of integrals irs = [] for itg_data_index, itg_data in enumerate(form_data.integral_data): @@ -221,15 +212,14 @@ def _compute_integral_ir( expression_ir = {} # Compute representation - entity_type = _entity_types[itg_data.integral_type] + integral_type = convert_to_integral_type(itg_data.integral_type) ufl_cell = itg_data.domain.ufl_cell() cell_type = basix_cell_from_string(ufl_cell.cellname()) tdim = ufl_cell.topological_dimension() assert all(tdim == itg.ufl_domain().topological_dimension() for itg in itg_data.integrals) expression_ir = { - "integral_type": itg_data.integral_type, - "entity_type": entity_type, + "integral_type": integral_type, "shape": (), "coordinate_element_hash": itg_data.domain.ufl_coordinate_element().basix_hash(), } @@ -251,12 +241,12 @@ def _compute_integral_ir( ] # Compute shape of element tensor - if expression_ir["integral_type"] == "interior_facet": - expression_ir["tensor_shape"] = [2 * dim for dim in argument_dimensions] - else: - expression_ir["tensor_shape"] = argument_dimensions - - integral_type = itg_data.integral_type + match expression_ir["integral_type"]: + case IntegralType(codim=1, num_neighbours=2): + # If interior facet integral, double the shape in the first index + expression_ir["tensor_shape"] = [2 * dim for dim in argument_dimensions] + case _: + expression_ir["tensor_shape"] = argument_dimensions # Group integrands with the same quadrature rule grouped_integrands: dict[ @@ -282,14 +272,15 @@ def _compute_integral_ir( # prescribed in certain cases. degree = md["quadrature_degree"] - if "facet" in integral_type: - facet_types = basix.cell.subentity_types(cell_type)[-2] - assert len(set(facet_types)) == 1 - cell_type = facet_types[0] - elif integral_type == "ridge": - ridge_types = basix.cell.subentity_types(cell_type)[-3] - assert len(set(ridge_types)) == 1 - cell_type = ridge_types[0] + match integral_type: + case IntegralType(codim=1): + facet_types = basix.cell.subentity_types(cell_type)[-2] + assert len(set(facet_types)) == 1 + cell_type = facet_types[0] + case IntegralType(codim=2): + ridge_types = basix.cell.subentity_types(cell_type)[-3] + assert len(set(ridge_types)) == 1 + cell_type = ridge_types[0] if degree > 1: warnings.warn( @@ -360,7 +351,13 @@ def _compute_integral_ir( index_to_coeff = sorted([(v, k) for k, v in coefficient_numbering.items()]) offsets = {} - width = 2 if integral_type in ("interior_facet") else 1 + match integral_type: + case IntegralType(codim=1, num_neighbours=2): + # Double width for interior facet integrals + width = 2 + case _: + width = 1 + _offset = 0 for k, el in zip(index_to_coeff, form_data.coefficient_elements): offsets[k[1]] = _offset @@ -387,8 +384,7 @@ def _compute_integral_ir( # Build more specific intermediate representation integral_ir = compute_integral_ir( itg_data.domain.ufl_cell(), - itg_data.integral_type, - expression_ir["entity_type"], + integral_type, integrand_map, expression_ir["tensor_shape"], options, @@ -454,14 +450,20 @@ def _compute_form_ir( # Store names of integrals and subdomain_ids for this form, grouped # by integral types since form points to all integrals it contains, # it has to know their names for codegen phase - ufcx_integral_types = ("cell", "exterior_facet", "interior_facet", "vertex", "ridge") + ufcx_integral_types = ( + IntegralType(codim=0), + IntegralType(codim=1, num_neighbours=1), + IntegralType(codim=1, num_neighbours=2), + IntegralType(codim=2), + IntegralType(codim=-1), + ) ir["subdomain_ids"] = {itg_type: [] for itg_type in ufcx_integral_types} ir["integral_names"] = {itg_type: [] for itg_type in ufcx_integral_types} ir["integral_domains"] = {itg_type: [] for itg_type in ufcx_integral_types} for itg_index, itg_data in enumerate(form_data.integral_data): # UFL is using "otherwise" for default integrals (over whole mesh) # but FFCx needs integers, so otherwise = -1 - integral_type = itg_data.integral_type + integral_type = convert_to_integral_type(itg_data.integral_type) subdomain_ids = [sid if sid != "otherwise" else -1 for sid in itg_data.subdomain_id] if min(subdomain_ids) < -1: @@ -563,12 +565,11 @@ def _compute_expression_ir( # Copy offsets also into IR base_ir["coefficient_offsets"] = offsets - base_ir["integral_type"] = "expression" if cell is not None: if (tdim := cell.topological_dimension()) == (pdim := points.shape[1]): - base_ir["entity_type"] = "cell" + base_ir["integral_type"] = IntegralType(codim=0, is_expression=True) elif tdim - 1 == pdim: - base_ir["entity_type"] = "facet" + base_ir["integral_type"] = IntegralType(codim=1, is_expression=True) else: raise ValueError( f"Expression on domain with topological dimension {tdim}" @@ -576,7 +577,7 @@ def _compute_expression_ir( ) else: # For spatially invariant expressions, all expressions are evaluated in the cell - base_ir["entity_type"] = "cell" + base_ir["integral_type"] = IntegralType(codim=0, is_expression=True) # Build offsets for Constants original_constant_offsets = {} @@ -605,7 +606,6 @@ def _compute_expression_ir( expression_ir = compute_integral_ir( cell, base_ir["integral_type"], - base_ir["entity_type"], integrands, tensor_shape, options, diff --git a/ffcx/ir/representationutils.py b/ffcx/ir/representationutils.py index 54c731db9..3201c277b 100644 --- a/ffcx/ir/representationutils.py +++ b/ffcx/ir/representationutils.py @@ -1,4 +1,4 @@ -# Copyright (C) 2012-2017 Marie Rognes +# Copyright (C) 2012-2025 Marie Rognes and Jørgen S. Dokken # # This file is part of FFCx.(https://www.fenicsproject.org) # @@ -8,10 +8,13 @@ import hashlib import itertools import logging +import typing -import numpy as np import ufl +import numpy as np +from basix.ufl import _ElementBase +from ffcx.definitions import IntegralType from ffcx.element_interface import ( create_quadrature, map_edge_points, @@ -55,77 +58,80 @@ def id(self): def create_quadrature_points_and_weights( - integral_type, cell, degree, rule, elements, use_tensor_product=False + integral_type: IntegralType, + cell: ufl.Cell, + degree: int, + rule: str, + elements: list[_ElementBase], + use_tensor_product: bool = False, ): """Create quadrature rule and return points and weights.""" pts = {} wts = {} tensor_factors = {} - if integral_type == "cell": - cell_name = cell.cellname() - if cell_name in ["quadrilateral", "hexahedron"] and use_tensor_product: - if cell_name == "quadrilateral": - tensor_factors[cell_name] = [ - create_quadrature("interval", degree, rule, elements) for _ in range(2) - ] - elif cell_name == "hexahedron": - tensor_factors[cell_name] = [ - create_quadrature("interval", degree, rule, elements) for _ in range(3) - ] - pts[cell_name] = np.array( - [ - tuple(i[0] for i in p) - for p in itertools.product(*[f[0] for f in tensor_factors[cell_name]]) - ] - ) - wts[cell_name] = np.array( - [np.prod(p) for p in itertools.product(*[f[1] for f in tensor_factors[cell_name]])] - ) - else: - pts[cell_name], wts[cell_name] = create_quadrature(cell_name, degree, rule, elements) - elif integral_type in ufl.measure.facet_integral_types: - for ft in cell.facet_types(): - pts[ft.cellname()], wts[ft.cellname()] = create_quadrature( - ft.cellname(), - degree, - rule, - elements, + match integral_type: + case IntegralType(is_expression=True): + pass + case IntegralType(codim=0): + cell_name = cell.cellname() + if cell_name in ["quadrilateral", "hexahedron"] and use_tensor_product: + if cell_name == "quadrilateral": + tensor_factors[cell_name] = [ + create_quadrature("interval", degree, rule, elements) for _ in range(2) + ] + elif cell_name == "hexahedron": + tensor_factors[cell_name] = [ + create_quadrature("interval", degree, rule, elements) for _ in range(3) + ] + pts[cell_name] = np.array( + [ + tuple(i[0] for i in p) # type: ignore + for p in itertools.product(*[f[0] for f in tensor_factors[cell_name]]) # type: ignore + ] + ) + wts[cell_name] = np.array( + [ + np.prod(p) # type: ignore + for p in itertools.product(*[f[1] for f in tensor_factors[cell_name]]) # type: ignore + ] + ) + else: + pts[cell_name], wts[cell_name] = create_quadrature( # type: ignore + cell_name, degree, rule, elements + ) + case IntegralType(codim=1): + for ft in cell.facet_types(): + pts[ft.cellname()], wts[ft.cellname()] = create_quadrature( # type: ignore + ft.cellname(), + degree, + rule, + elements, + ) + case IntegralType(codim=2): + for rt in cell.ridge_types(): + pts[rt.cellname()], wts[rt.cellname()] = create_quadrature( # type: ignore + rt.cellname(), + degree, + rule, + elements, + ) + case IntegralType(codim=-1): + pts["vertex"], wts["vertex"] = create_quadrature("vertex", degree, rule, elements) # type: ignore + case _: + raise RuntimeError( + f"Unknown integral_type {integral_type} in quadrature rule creation." ) - elif integral_type in ufl.measure.ridge_integral_types: - for rt in cell.ridge_types(): - pts[rt.cellname()], wts[rt.cellname()] = create_quadrature( - rt.cellname(), - degree, - rule, - elements, - ) - elif integral_type in ufl.measure.point_integral_types: - pts["vertex"], wts["vertex"] = create_quadrature("vertex", degree, rule, elements) - elif integral_type == "expression": - pass - else: - logger.exception(f"Unknown integral type: {integral_type}") - return pts, wts, tensor_factors -def integral_type_to_entity_dim(integral_type, tdim): +def integral_type_to_entity_dim(integral_type: IntegralType, tdim: int): """Given integral_type and domain tdim, return the tdim of the integration entity.""" - if integral_type == "cell": - entity_dim = tdim - elif integral_type in ufl.measure.facet_integral_types: - entity_dim = tdim - 1 - elif integral_type in ufl.measure.ridge_integral_types: - entity_dim = tdim - 2 - elif integral_type in ufl.measure.point_integral_types: - entity_dim = 0 - elif integral_type in ufl.custom_integral_types: - entity_dim = tdim - elif integral_type == "expression": - entity_dim = tdim + if integral_type.codim == -1: + return 0 + elif integral_type.codim >= 0: + return tdim - integral_type.codim else: - raise RuntimeError(f"Unknown integral_type: {integral_type}") - return entity_dim + raise NotImplementedError(f"Cannot handle codim {integral_type.codim} < -1") def map_integral_points(points, integral_type, cell, entity): diff --git a/test/test_add_mode.py b/test/test_add_mode.py index f57fde78b..e205e26b8 100644 --- a/test/test_add_mode.py +++ b/test/test_add_mode.py @@ -56,7 +56,7 @@ def test_additive_facet_integral(dtype, compile_args): form0 = compiled_forms[0] integral_offsets = form0.form_integral_offsets - ex = module.lib.exterior_facet + ex = module.lib.facet assert integral_offsets[ex + 1] - integral_offsets[ex] == 1 integral_id = form0.form_integral_ids[integral_offsets[ex]] assert integral_id == -1 diff --git a/test/test_jit_forms.py b/test/test_jit_forms.py index 5d01571c4..0a29956d4 100644 --- a/test/test_jit_forms.py +++ b/test/test_jit_forms.py @@ -391,7 +391,7 @@ def test_subdomains(compile_args): form3 = compiled_forms[3] offsets = form3.form_integral_offsets assert offsets[cell + 1] - offsets[cell] == 0 - exf = module.lib.exterior_facet + exf = module.lib.facet ids = [form3.form_integral_ids[j] for j in range(offsets[exf], offsets[exf + 1])] assert ids[0] == 0 and ids[1] == 210 @@ -1028,7 +1028,7 @@ def test_facet_vertex_quadrature(compile_args): solutions = [] for form in compiled_forms: offsets = form.form_integral_offsets - exf = module.lib.exterior_facet + exf = module.lib.facet assert offsets[exf + 1] - offsets[exf] == 1 default_integral = form.form_integrals[offsets[exf]] @@ -1311,7 +1311,7 @@ def test_ds_prism(compile_args, dtype): offsets = form0.form_integral_offsets cell = module.lib.cell - exterior_facet = module.lib.exterior_facet + exterior_facet = module.lib.facet interior_facet = module.lib.interior_facet assert offsets[cell + 1] - offsets[cell] == 0 assert offsets[exterior_facet + 1] - offsets[exterior_facet] == 2 From 1aca402878e5cb3fca14b81651801f48c718413c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 09:59:12 +0000 Subject: [PATCH 02/10] Typing on mac differs from my typing locally --- ffcx/analysis.py | 6 +++--- ffcx/codegeneration/access.py | 4 ++-- ffcx/codegeneration/definitions.py | 4 ++-- ffcx/ir/elementtables.py | 2 +- ffcx/ir/representation.py | 6 +++--- 5 files changed, 11 insertions(+), 11 deletions(-) diff --git a/ffcx/analysis.py b/ffcx/analysis.py index dba9c6be3..31d214012 100644 --- a/ffcx/analysis.py +++ b/ffcx/analysis.py @@ -93,8 +93,8 @@ def analyze_ufl_objects( form_data = tuple(_analyze_form(form, scalar_type) for form in forms) for data in form_data: - elements += data.unique_sub_elements - coordinate_elements += data.coordinate_elements + elements += data.unique_sub_elements # type: ignore + coordinate_elements += data.coordinate_elements # type: ignore for original_expression, points in expressions: elements += ufl.algorithms.extract_elements(original_expression) @@ -190,7 +190,7 @@ def _analyze_form( # Determine unique quadrature degree and quadrature scheme # per each integral data - for id, integral_data in enumerate(form_data.integral_data): + for id, integral_data in enumerate(form_data.integral_data): # type: ignore # Iterate through groups of integral data. There is one integral # data for all integrals with same domain, itype, subdomain_id # (but possibly different metadata). diff --git a/ffcx/codegeneration/access.py b/ffcx/codegeneration/access.py index 3125a71ef..a70d1478f 100644 --- a/ffcx/codegeneration/access.py +++ b/ffcx/codegeneration/access.py @@ -64,7 +64,7 @@ def get( """Format a terminal.""" e = mt.terminal # Call appropriate handler, depending on the type of e - handler = self.call_lookup.get(type(e), False) + handler = self.call_lookup.get(type(e), False) # type: ignore if not handler: # Look for parent class types instead @@ -73,7 +73,7 @@ def get( handler = self.call_lookup[k] break if handler: - return handler(mt, tabledata, quadrature_rule) # type: ignore + return handler(mt, tabledata, quadrature_rule) else: raise RuntimeError(f"Not handled: {type(e)}") diff --git a/ffcx/codegeneration/definitions.py b/ffcx/codegeneration/definitions.py index 3f4b13220..7523a2e64 100644 --- a/ffcx/codegeneration/definitions.py +++ b/ffcx/codegeneration/definitions.py @@ -101,13 +101,13 @@ def get( ttype = ttype.__bases__[0] # Get the handler from the lookup, or None if not found - handler = self.handler_lookup.get(ttype) + handler = self.handler_lookup.get(ttype) # type: ignore if handler is None: raise NotImplementedError(f"No handler for terminal type: {ttype}") # Call the handler - return handler(mt, tabledata, quadrature_rule, access) + return handler(mt, tabledata, quadrature_rule, access) # type: ignore def coefficient( self, diff --git a/ffcx/ir/elementtables.py b/ffcx/ir/elementtables.py index e6582f8c9..35e0f367d 100644 --- a/ffcx/ir/elementtables.py +++ b/ffcx/ir/elementtables.py @@ -230,7 +230,7 @@ def get_modified_terminal_element(mt) -> ModifiedTerminalElement | None: raise RuntimeError("Global derivatives of reference values not defined.") elif ld and not mt.reference_value: raise RuntimeError("Local derivatives of global values not defined.") - element = mt.terminal.ufl_function_space().ufl_element() + element = mt.terminal.ufl_function_space().ufl_element() # type: ignore fc = mt.flat_component elif isinstance(mt.terminal, ufl.classes.SpatialCoordinate): if mt.reference_value: diff --git a/ffcx/ir/representation.py b/ffcx/ir/representation.py index f5fe77d82..d0f49affa 100644 --- a/ffcx/ir/representation.py +++ b/ffcx/ir/representation.py @@ -137,10 +137,10 @@ def compute_ir( integral_names = {} form_names = {} for fd_index, fd in enumerate(analysis.form_data): - form_names[fd_index] = naming.form_name(fd.original_form, fd_index, prefix) - for itg_index, itg_data in enumerate(fd.integral_data): + form_names[fd_index] = naming.form_name(fd.original_form, fd_index, prefix) # type: ignore + for itg_index, itg_data in enumerate(fd.integral_data): # type: ignore integral_names[(fd_index, itg_index)] = naming.integral_name( - fd.original_form, + fd.original_form, # type: ignore itg_data.integral_type, fd_index, itg_data.subdomain_id, From 97fc82d2bcda88657a162778bb2f973726a199a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 10:34:05 +0000 Subject: [PATCH 03/10] Fix expr typing --- ffcx/codegeneration/jit.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ffcx/codegeneration/jit.py b/ffcx/codegeneration/jit.py index 6eb5dbb8f..22b01ad9b 100644 --- a/ffcx/codegeneration/jit.py +++ b/ffcx/codegeneration/jit.py @@ -229,7 +229,7 @@ def compile_forms( def compile_expressions( - expressions: list[tuple[ufl.Expr, npt.NDArray[np.floating]]], + expressions: list[tuple[ufl.core.expr.Expr, npt.NDArray[np.floating]]], options: dict = {}, cache_dir: Path | None = None, timeout: int = 10, From 82fac46123e8bd90a006987d92e7a3ec00ab577e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 10:42:16 +0000 Subject: [PATCH 04/10] Ruf fformatting --- ffcx/codegeneration/definitions.py | 2 +- ffcx/ir/elementtables.py | 7 ++++++- ffcx/ir/representationutils.py | 11 +++++------ 3 files changed, 12 insertions(+), 8 deletions(-) diff --git a/ffcx/codegeneration/definitions.py b/ffcx/codegeneration/definitions.py index 7523a2e64..32eca9bc5 100644 --- a/ffcx/codegeneration/definitions.py +++ b/ffcx/codegeneration/definitions.py @@ -107,7 +107,7 @@ def get( raise NotImplementedError(f"No handler for terminal type: {ttype}") # Call the handler - return handler(mt, tabledata, quadrature_rule, access) # type: ignore + return handler(mt, tabledata, quadrature_rule, access) # type: ignore def coefficient( self, diff --git a/ffcx/ir/elementtables.py b/ffcx/ir/elementtables.py index 35e0f367d..57b2dc25a 100644 --- a/ffcx/ir/elementtables.py +++ b/ffcx/ir/elementtables.py @@ -382,7 +382,12 @@ def build_optimized_tables( element_number = element_numbers[element] name = generate_psi_table_name( - quadrature_rule, element_number, avg, integral_type, local_derivatives, flat_component # type: ignore + quadrature_rule, + element_number, + avg, + integral_type, + local_derivatives, + flat_component, # type: ignore ) # FIXME - currently just recalculate the tables every time, diff --git a/ffcx/ir/representationutils.py b/ffcx/ir/representationutils.py index 3201c277b..07b2f186a 100644 --- a/ffcx/ir/representationutils.py +++ b/ffcx/ir/representationutils.py @@ -8,10 +8,9 @@ import hashlib import itertools import logging -import typing -import ufl import numpy as np +import ufl from basix.ufl import _ElementBase from ffcx.definitions import IntegralType @@ -96,12 +95,12 @@ def create_quadrature_points_and_weights( ] ) else: - pts[cell_name], wts[cell_name] = create_quadrature( # type: ignore + pts[cell_name], wts[cell_name] = create_quadrature( # type: ignore cell_name, degree, rule, elements ) case IntegralType(codim=1): for ft in cell.facet_types(): - pts[ft.cellname()], wts[ft.cellname()] = create_quadrature( # type: ignore + pts[ft.cellname()], wts[ft.cellname()] = create_quadrature( # type: ignore ft.cellname(), degree, rule, @@ -109,14 +108,14 @@ def create_quadrature_points_and_weights( ) case IntegralType(codim=2): for rt in cell.ridge_types(): - pts[rt.cellname()], wts[rt.cellname()] = create_quadrature( # type: ignore + pts[rt.cellname()], wts[rt.cellname()] = create_quadrature( # type: ignore rt.cellname(), degree, rule, elements, ) case IntegralType(codim=-1): - pts["vertex"], wts["vertex"] = create_quadrature("vertex", degree, rule, elements) # type: ignore + pts["vertex"], wts["vertex"] = create_quadrature("vertex", degree, rule, elements) # type: ignore case _: raise RuntimeError( f"Unknown integral_type {integral_type} in quadrature rule creation." From 1db1251e76400f5a43293a304c4162fb94c0165d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 10:44:29 +0000 Subject: [PATCH 05/10] Move ignore after formatting --- ffcx/ir/elementtables.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ffcx/ir/elementtables.py b/ffcx/ir/elementtables.py index 57b2dc25a..11a2b3e17 100644 --- a/ffcx/ir/elementtables.py +++ b/ffcx/ir/elementtables.py @@ -384,10 +384,10 @@ def build_optimized_tables( name = generate_psi_table_name( quadrature_rule, element_number, - avg, + avg, # type: ignore integral_type, local_derivatives, - flat_component, # type: ignore + flat_component, ) # FIXME - currently just recalculate the tables every time, From 52eb6256fb94d0cd884e7a582fffabe4717b58ae Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Wed, 10 Sep 2025 12:49:31 +0200 Subject: [PATCH 06/10] Swap order to mimic old code --- ffcx/ir/representation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ffcx/ir/representation.py b/ffcx/ir/representation.py index d0f49affa..c26d0af46 100644 --- a/ffcx/ir/representation.py +++ b/ffcx/ir/representation.py @@ -454,8 +454,8 @@ def _compute_form_ir( IntegralType(codim=0), IntegralType(codim=1, num_neighbours=1), IntegralType(codim=1, num_neighbours=2), - IntegralType(codim=2), IntegralType(codim=-1), + IntegralType(codim=2), ) ir["subdomain_ids"] = {itg_type: [] for itg_type in ufcx_integral_types} ir["integral_names"] = {itg_type: [] for itg_type in ufcx_integral_types} From b57309e9aabf204bb07de298294710ffbc1f7972 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 10:51:59 +0000 Subject: [PATCH 07/10] Fix some more access call --- ffcx/codegeneration/access.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ffcx/codegeneration/access.py b/ffcx/codegeneration/access.py index a70d1478f..fce221566 100644 --- a/ffcx/codegeneration/access.py +++ b/ffcx/codegeneration/access.py @@ -303,7 +303,7 @@ def facet_orientation(self, mt, tabledata, num_points): raise RuntimeError(f"Unhandled cell types {cellname}.") table = L.Symbol(f"{cellname}_facet_orientation", dtype=L.DataType.INT) - facet = self.symbols.entity("facet", mt.restriction) + facet = self.symbols.entity(IntegralType(codim=1), mt.restriction) return table[facet] def cell_vertices(self, mt, tabledata, num_points): @@ -398,7 +398,7 @@ def facet_edge_vectors(self, mt, tabledata, num_points): num_scalar_dofs = scalar_element.dim # Get edge vertices - facet = self.symbols.entity("facet", mt.restriction) + facet = self.symbols.entity(IntegralType(codim=1), mt.restriction) facet_edge = mt.component[0] facet_edge_vertices = L.Symbol(f"{cellname}_facet_edge_vertices", dtype=L.DataType.INT) vertex0 = facet_edge_vertices[facet][facet_edge][0] From b77d383793ff508f4e86bb94c8f23b362abac36e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 10:56:51 +0000 Subject: [PATCH 08/10] Change dolfinx branch on FFCx --- .github/workflows/dolfinx-tests.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/dolfinx-tests.yml b/.github/workflows/dolfinx-tests.yml index 41c1b4586..37462b54f 100644 --- a/.github/workflows/dolfinx-tests.yml +++ b/.github/workflows/dolfinx-tests.yml @@ -55,7 +55,7 @@ jobs: with: path: ./dolfinx repository: FEniCS/dolfinx - ref: main + ref: dokken/itg_type - name: Get DOLFINx source (specified branch/tag) if: github.event_name == 'workflow_dispatch' uses: actions/checkout@v5 From 8128617c857a55dc40bf6f6d182d35758557e9b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Wed, 10 Sep 2025 12:59:01 +0200 Subject: [PATCH 09/10] Add back comments --- ffcx/ir/elementtables.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ffcx/ir/elementtables.py b/ffcx/ir/elementtables.py index 11a2b3e17..bfd9070e2 100644 --- a/ffcx/ir/elementtables.py +++ b/ffcx/ir/elementtables.py @@ -98,10 +98,13 @@ def get_ffcx_table_values( match integral_type: case IntegralType(is_custom=True): + # Use quadrature points on cell for analysis in custom integral types assert integral_type.codim == 0 integral_type = IntegralType(codim=0) assert not avg case IntegralType(is_expression=True): + # FFCx tables for expression are generated as either interior cell points + # or points on a facet assert integral_type.codim in (0, 1) assert integral_type.num_neighbours == 1 From e2d03ae376697cc3a52552259dc6d8dbf77a0af8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Wed, 10 Sep 2025 13:11:02 +0000 Subject: [PATCH 10/10] Rename num_neighbours to num_cells --- ffcx/codegeneration/C/form.py | 4 ++-- ffcx/definitions.py | 14 +++++++------- ffcx/ir/elementtables.py | 4 ++-- ffcx/ir/representation.py | 8 ++++---- 4 files changed, 15 insertions(+), 15 deletions(-) diff --git a/ffcx/codegeneration/C/form.py b/ffcx/codegeneration/C/form.py index 8d33a4044..242c7f49d 100644 --- a/ffcx/codegeneration/C/form.py +++ b/ffcx/codegeneration/C/form.py @@ -122,8 +122,8 @@ def generator(ir: FormIR, options): # Note: the order of this list is defined by the enum ufcx_integral_type in ufcx.h for itg_type in ( IntegralType(codim=0), - IntegralType(codim=1, num_neighbours=1), - IntegralType(codim=1, num_neighbours=2), + IntegralType(codim=1, num_cells=1), + IntegralType(codim=1, num_cells=2), IntegralType(codim=-1), IntegralType(codim=2), ): diff --git a/ffcx/definitions.py b/ffcx/definitions.py index 7c304e971..d757049e3 100644 --- a/ffcx/definitions.py +++ b/ffcx/definitions.py @@ -13,7 +13,7 @@ class IntegralType(NamedTuple): """Class for storing information about an integral type.""" codim: int # Codimension of integral with respect to the topological dimension of the domain - num_neighbours: int = 1 # Number of neighbouring cells + num_cells: int = 1 # Number of cells the integral is over is_expression: bool = ( False # True if the integral is an Expression (integrand) rather than an integral ) @@ -24,16 +24,16 @@ def convert_to_integral_type(integral_type: str) -> IntegralType: """Convert UFL integral type string to IntegralType.""" match integral_type: case "cell": - return IntegralType(codim=0, num_neighbours=1) + return IntegralType(codim=0, num_cells=1) case "exterior_facet": - return IntegralType(codim=1, num_neighbours=1) + return IntegralType(codim=1, num_cells=1) case "interior_facet": - return IntegralType(codim=1, num_neighbours=2) + return IntegralType(codim=1, num_cells=2) case "ridge": - return IntegralType(codim=2, num_neighbours=1) + return IntegralType(codim=2, num_cells=1) case "vertex": - return IntegralType(codim=-1, num_neighbours=1) + return IntegralType(codim=-1, num_cells=1) case "custom": - return IntegralType(codim=0, num_neighbours=1, is_custom=True) + return IntegralType(codim=0, num_cells=1, is_custom=True) case _: raise ValueError(f"Unknown integral type:{integral_type}") diff --git a/ffcx/ir/elementtables.py b/ffcx/ir/elementtables.py index bfd9070e2..26ad28769 100644 --- a/ffcx/ir/elementtables.py +++ b/ffcx/ir/elementtables.py @@ -106,7 +106,7 @@ def get_ffcx_table_values( # FFCx tables for expression are generated as either interior cell points # or points on a facet assert integral_type.codim in (0, 1) - assert integral_type.num_neighbours == 1 + assert integral_type.num_cells == 1 if avg in ("cell", "facet"): # Redefine points to compute average tables @@ -409,7 +409,7 @@ def build_optimized_tables( # needed because a cell may see its sub-entities as being oriented # differently to their global orientation if ( - (integral_type.num_neighbours > 1 and integral_type.codim == 1) + (integral_type.num_cells > 1 and integral_type.codim == 1) or integral_type.codim == 2 or (is_mixed_dim and codim == 0) ): diff --git a/ffcx/ir/representation.py b/ffcx/ir/representation.py index c26d0af46..6cf5bf586 100644 --- a/ffcx/ir/representation.py +++ b/ffcx/ir/representation.py @@ -242,7 +242,7 @@ def _compute_integral_ir( # Compute shape of element tensor match expression_ir["integral_type"]: - case IntegralType(codim=1, num_neighbours=2): + case IntegralType(codim=1, num_cells=2): # If interior facet integral, double the shape in the first index expression_ir["tensor_shape"] = [2 * dim for dim in argument_dimensions] case _: @@ -352,7 +352,7 @@ def _compute_integral_ir( index_to_coeff = sorted([(v, k) for k, v in coefficient_numbering.items()]) offsets = {} match integral_type: - case IntegralType(codim=1, num_neighbours=2): + case IntegralType(codim=1, num_cells=2): # Double width for interior facet integrals width = 2 case _: @@ -452,8 +452,8 @@ def _compute_form_ir( # it has to know their names for codegen phase ufcx_integral_types = ( IntegralType(codim=0), - IntegralType(codim=1, num_neighbours=1), - IntegralType(codim=1, num_neighbours=2), + IntegralType(codim=1, num_cells=1), + IntegralType(codim=1, num_cells=2), IntegralType(codim=-1), IntegralType(codim=2), )