From e5e8209da7b76368417fe8ef78f4173d9ed2c2ed Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Wed, 19 Aug 2026 08:19:54 +0200 Subject: [PATCH 01/11] Add scalar geometry option --- ffcx_backends/cpp.py | 29 ++++++++++++++++++++--------- 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/ffcx_backends/cpp.py b/ffcx_backends/cpp.py index 6283ddc..40877b9 100644 --- a/ffcx_backends/cpp.py +++ b/ffcx_backends/cpp.py @@ -19,12 +19,12 @@ logger = logging.getLogger("ffcx") -def dtype_to_cpp_type(dtype: L.DataType, scalar_type: str, real_type: str) -> str: +def dtype_to_cpp_type(dtype: L.DataType, scalar_type: str, geometry_type: str) -> str: """Map L.DataType to C++ type.""" if dtype == L.DataType.SCALAR: return scalar_type elif dtype == L.DataType.REAL: - return real_type + return geometry_type elif dtype == L.DataType.INT: return "std::int32_t" elif dtype == L.DataType.BOOL: @@ -33,6 +33,11 @@ def dtype_to_cpp_type(dtype: L.DataType, scalar_type: str, real_type: str) -> st raise ValueError(f"Invalid datatype: {dtype}") +def geometry_type_name(options: dict[str, Any]) -> str: + """C++ template parameter the geometry-derived temporaries are emitted in.""" + return "T" if options.get("scalar_geometry", False) else "U" + + class Formatter: """Format FFCx nodes into C++.""" @@ -76,10 +81,16 @@ def build_initializer_lists(values: npt.NDArray) -> str: arr += "}" return arr - def __init__(self, scalar: Any) -> None: - """Initialise.""" + def __init__(self, scalar: Any, geometry_type: str = "U") -> None: + """Initialise. + + Args: + scalar: Scalar type of the element tensor. + geometry_type: C++ template parameter the geometry-derived temporaries are + emitted in, see :func:`geometry_type_name`. + """ self.scalar_type = "T" - self.real_type = "U" + self.geometry_type = geometry_type @functools.singledispatchmethod def __call__(self, obj: L.LNode) -> str: @@ -122,7 +133,7 @@ def format_array_decl(self, arr: L.ArrayDecl) -> str: dtype = arr.symbol.dtype assert dtype is not None - typename = dtype_to_cpp_type(dtype, self.scalar_type, self.real_type) + typename = dtype_to_cpp_type(dtype, self.scalar_type, self.geometry_type) symbol = self(arr.symbol) dims = "".join([f"[{i}]" for i in arr.sizes]) @@ -152,7 +163,7 @@ def format_variable_decl(self, v: L.VariableDecl) -> str: val = self(v.value) symbol = self(v.symbol) assert v.symbol.dtype - typename = dtype_to_cpp_type(v.symbol.dtype, self.scalar_type, self.real_type) + typename = dtype_to_cpp_type(v.symbol.dtype, self.scalar_type, self.geometry_type) return f"{typename} {symbol} = {val};\n" @__call__.register @@ -329,7 +340,7 @@ def generator(ir: ExpressionIR, options: dict[str, int | float | npt.DTypeLike]) d["factory_name"] = factory_name parts = eg.generate() - formatter = Formatter(options["scalar_type"]) + formatter = Formatter(options["scalar_type"], geometry_type_name(options)) d["tabulate_expression"] = formatter(parts) if len(ir.original_coefficient_positions) > 0: @@ -442,7 +453,7 @@ def generator( parts = ig.generate(domain) # Format code as string - formatter = Formatter(options["scalar_type"]) + formatter = Formatter(options["scalar_type"], geometry_type_name(options)) body = formatter(parts) # Generate generic FFCx code snippets and add specific parts From 3a5786234502ee9867719bf4d48172bab4d9bba4 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 25 Aug 2026 16:44:47 +0200 Subject: [PATCH 02/11] Update --- ffcx_backends/cpp.py | 54 +++++++++++++++++-------------------- test/test_cpp.py | 63 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 87 insertions(+), 30 deletions(-) diff --git a/ffcx_backends/cpp.py b/ffcx_backends/cpp.py index 40877b9..d2c2fce 100644 --- a/ffcx_backends/cpp.py +++ b/ffcx_backends/cpp.py @@ -4,7 +4,7 @@ import logging import pprint import textwrap -from typing import Any, ClassVar +from typing import ClassVar import basix import ffcx.codegeneration.lnodes as L # noqa @@ -19,25 +19,6 @@ logger = logging.getLogger("ffcx") -def dtype_to_cpp_type(dtype: L.DataType, scalar_type: str, geometry_type: str) -> str: - """Map L.DataType to C++ type.""" - if dtype == L.DataType.SCALAR: - return scalar_type - elif dtype == L.DataType.REAL: - return geometry_type - elif dtype == L.DataType.INT: - return "std::int32_t" - elif dtype == L.DataType.BOOL: - return "bool" - else: - raise ValueError(f"Invalid datatype: {dtype}") - - -def geometry_type_name(options: dict[str, Any]) -> str: - """C++ template parameter the geometry-derived temporaries are emitted in.""" - return "T" if options.get("scalar_geometry", False) else "U" - - class Formatter: """Format FFCx nodes into C++.""" @@ -81,16 +62,29 @@ def build_initializer_lists(values: npt.NDArray) -> str: arr += "}" return arr - def __init__(self, scalar: Any, geometry_type: str = "U") -> None: + def __init__(self, scalar_geometry: bool = False) -> None: """Initialise. Args: - scalar: Scalar type of the element tensor. - geometry_type: C++ template parameter the geometry-derived temporaries are - emitted in, see :func:`geometry_type_name`. + scalar_geometry: Force scalar type == geometry type in the kernel body. """ - self.scalar_type = "T" - self.geometry_type = geometry_type + self._scalar_type = "T" + self._geometry_type = "U" + if scalar_geometry: + self._geometry_type = self._scalar_type + + def dtype_to_cpp_type(self, dtype: L.DataType) -> str: + """Map L.DataType to C++ type.""" + if dtype == L.DataType.SCALAR: + return self._scalar_type + elif dtype == L.DataType.REAL: + return self._geometry_type + elif dtype == L.DataType.INT: + return "std::int32_t" + elif dtype == L.DataType.BOOL: + return "bool" + else: + raise ValueError(f"Invalid datatype: {dtype}") @functools.singledispatchmethod def __call__(self, obj: L.LNode) -> str: @@ -133,7 +127,7 @@ def format_array_decl(self, arr: L.ArrayDecl) -> str: dtype = arr.symbol.dtype assert dtype is not None - typename = dtype_to_cpp_type(dtype, self.scalar_type, self.geometry_type) + typename = self.dtype_to_cpp_type(dtype) symbol = self(arr.symbol) dims = "".join([f"[{i}]" for i in arr.sizes]) @@ -163,7 +157,7 @@ def format_variable_decl(self, v: L.VariableDecl) -> str: val = self(v.value) symbol = self(v.symbol) assert v.symbol.dtype - typename = dtype_to_cpp_type(v.symbol.dtype, self.scalar_type, self.geometry_type) + typename = self.dtype_to_cpp_type(v.symbol.dtype) return f"{typename} {symbol} = {val};\n" @__call__.register @@ -340,7 +334,7 @@ def generator(ir: ExpressionIR, options: dict[str, int | float | npt.DTypeLike]) d["factory_name"] = factory_name parts = eg.generate() - formatter = Formatter(options["scalar_type"], geometry_type_name(options)) + formatter = Formatter(options.get("scalar_geometry", False)) d["tabulate_expression"] = formatter(parts) if len(ir.original_coefficient_positions) > 0: @@ -453,7 +447,7 @@ def generator( parts = ig.generate(domain) # Format code as string - formatter = Formatter(options["scalar_type"], geometry_type_name(options)) + formatter = Formatter(options.get("scalar_geometry", False)) body = formatter(parts) # Generate generic FFCx code snippets and add specific parts diff --git a/test/test_cpp.py b/test/test_cpp.py index 57cd16c..d88e32b 100644 --- a/test/test_cpp.py +++ b/test/test_cpp.py @@ -1,8 +1,14 @@ +import re + import basix.ufl +import ffcx.codegeneration.lnodes as L # noqa: N812 +import pytest import ufl from ffcx.compiler import compile_ufl_objects from ffcx.options import get_options +from ffcx_backends.cpp import Formatter + def test_integral() -> None: element = basix.ufl.element("Lagrange", "triangle", 1) @@ -16,3 +22,60 @@ def test_integral() -> None: compiled_objects = compile_ufl_objects([a], opts) assert len(compiled_objects) == 2 + + +@pytest.mark.parametrize(("scalar_geometry", "expected"), [(False, "U"), (True, "T"), (None, "U")]) +def test_scalar_geometry(scalar_geometry: bool | None, expected: str) -> None: + """Geometry-derived temporaries are emitted in T only with ``scalar_geometry``. + + The kernel signature is unaffected: coordinate dofs are always read as ``U``, + the option only changes the type they are computed in. + """ + element = basix.ufl.element("Lagrange", "triangle", 1) + domain = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(2,))) + space = ufl.FunctionSpace(domain, element) + u, v = ufl.TrialFunction(space), ufl.TestFunction(space) + + a = (ufl.inner(u, v) + ufl.inner(ufl.grad(u), ufl.grad(v))) * ufl.dx + opts = get_options({"language": "ffcx_backends.cpp"}) + if scalar_geometry is not None: + opts["scalar_geometry"] = scalar_geometry + + code = compile_ufl_objects([a], opts)[0][0] + + # Quadrature weights, tabulated basis functions and the Jacobian are all + # geometry-derived, hence declared in the geometry type. + for declaration in [ + r"static const (\w+) weights_\w+\[", + r"static const (\w+) FE\w+\[", + # The Jacobian symbol carries a process-global counter, hence J. + r"^(\w+) J\d+_c0 = ", + ]: + matches = re.findall(declaration, code, flags=re.MULTILINE) + assert matches, f"no declaration matching {declaration} in generated code" + assert set(matches) == {expected} + + # The tabulate_tensor signature keeps the scalar/geometry split either way. + signature = re.search(r"static void tabulate_tensor\((.*?)\)\s*\{", code, flags=re.DOTALL) + assert signature is not None + assert "const U* RESTRICT coordinate_dofs" in signature.group(1) + assert "T* RESTRICT A" in signature.group(1) + + +@pytest.mark.parametrize(("scalar_geometry", "expected"), [(False, "U"), (True, "T")]) +def test_formatter_dtype_to_cpp_type(scalar_geometry: bool, expected: str) -> None: + """REAL follows the geometry type, the remaining datatypes do not.""" + formatter = Formatter(scalar_geometry) + + assert formatter.dtype_to_cpp_type(L.DataType.REAL) == expected + assert formatter.dtype_to_cpp_type(L.DataType.SCALAR) == "T" + assert formatter.dtype_to_cpp_type(L.DataType.INT) == "std::int32_t" + assert formatter.dtype_to_cpp_type(L.DataType.BOOL) == "bool" + + with pytest.raises(ValueError, match="Invalid datatype"): + formatter.dtype_to_cpp_type("not-a-datatype") # type: ignore[arg-type] + + +def test_formatter_default_geometry_type() -> None: + """Geometry stays in its own type unless asked otherwise.""" + assert Formatter().dtype_to_cpp_type(L.DataType.REAL) == "U" From 1199504c44f4193a9d157376f043d2f2e9159a83 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 25 Aug 2026 16:50:34 +0200 Subject: [PATCH 03/11] Remove static std::vector --- ffcx_backends/cpp.py | 2 +- test/test_cpp.py | 4 +++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/ffcx_backends/cpp.py b/ffcx_backends/cpp.py index d2c2fce..c7342a9 100644 --- a/ffcx_backends/cpp.py +++ b/ffcx_backends/cpp.py @@ -419,7 +419,7 @@ class {factory_name} }} // Data - static inline const std::vector enabled_coefficients{enabled_coefficients_init}; + static constexpr int enabled_coefficients[] = {enabled_coefficients_init}; static constexpr bool needs_facet_permutations = {needs_facet_permutations}; }}; diff --git a/test/test_cpp.py b/test/test_cpp.py index d88e32b..88bdc63 100644 --- a/test/test_cpp.py +++ b/test/test_cpp.py @@ -15,11 +15,13 @@ def test_integral() -> None: domain = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 1, shape=(2,))) space = ufl.FunctionSpace(domain, element) u, v = ufl.TrialFunction(space), ufl.TestFunction(space) + f = ufl.Coefficient(space) - a = (ufl.inner(u, v) + ufl.inner(ufl.grad(u), ufl.grad(v))) * ufl.dx + a = (ufl.inner(u, v) + f * ufl.inner(ufl.grad(u), ufl.grad(v))) * ufl.dx opts = get_options({"language": "ffcx_backends.cpp"}) compiled_objects = compile_ufl_objects([a], opts) + print(compiled_objects[0][0]) assert len(compiled_objects) == 2 From 7c6f6a8443cbca6949444c69239f5716f631249d Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 25 Aug 2026 16:53:16 +0200 Subject: [PATCH 04/11] Help poor mypy --- ffcx_backends/cpp.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/ffcx_backends/cpp.py b/ffcx_backends/cpp.py index d2c2fce..1631ae4 100644 --- a/ffcx_backends/cpp.py +++ b/ffcx_backends/cpp.py @@ -334,7 +334,8 @@ def generator(ir: ExpressionIR, options: dict[str, int | float | npt.DTypeLike]) d["factory_name"] = factory_name parts = eg.generate() - formatter = Formatter(options.get("scalar_geometry", False)) + scalar_geometry = True if options.get("scalar_geometry") else False + formatter = Formatter(scalar_geometry) d["tabulate_expression"] = formatter(parts) if len(ir.original_coefficient_positions) > 0: @@ -447,7 +448,8 @@ def generator( parts = ig.generate(domain) # Format code as string - formatter = Formatter(options.get("scalar_geometry", False)) + scalar_geometry = True if options.get("scalar_geometry") else False + formatter = Formatter(scalar_geometry) body = formatter(parts) # Generate generic FFCx code snippets and add specific parts From c4e7db5600ce1b4d15f4d5bb4e8061abde6a4297 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 25 Aug 2026 16:54:05 +0200 Subject: [PATCH 05/11] Remove debug print --- test/test_cpp.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_cpp.py b/test/test_cpp.py index 88bdc63..49d0e71 100644 --- a/test/test_cpp.py +++ b/test/test_cpp.py @@ -21,7 +21,6 @@ def test_integral() -> None: opts = get_options({"language": "ffcx_backends.cpp"}) compiled_objects = compile_ufl_objects([a], opts) - print(compiled_objects[0][0]) assert len(compiled_objects) == 2 From 76f57da53f0ed434e0ca6a1e89601e430f223a80 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 25 Aug 2026 16:58:11 +0200 Subject: [PATCH 06/11] simplify --- test/test_cpp.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/test/test_cpp.py b/test/test_cpp.py index d88e32b..2b7b99f 100644 --- a/test/test_cpp.py +++ b/test/test_cpp.py @@ -72,9 +72,6 @@ def test_formatter_dtype_to_cpp_type(scalar_geometry: bool, expected: str) -> No assert formatter.dtype_to_cpp_type(L.DataType.INT) == "std::int32_t" assert formatter.dtype_to_cpp_type(L.DataType.BOOL) == "bool" - with pytest.raises(ValueError, match="Invalid datatype"): - formatter.dtype_to_cpp_type("not-a-datatype") # type: ignore[arg-type] - def test_formatter_default_geometry_type() -> None: """Geometry stays in its own type unless asked otherwise.""" From 721cd8584bac535db8372252e4140924cacd8e19 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 25 Aug 2026 21:41:40 +0200 Subject: [PATCH 07/11] Use {0} --- ffcx_backends/cpp.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/ffcx_backends/cpp.py b/ffcx_backends/cpp.py index 8a1189d..59b18a8 100644 --- a/ffcx_backends/cpp.py +++ b/ffcx_backends/cpp.py @@ -458,7 +458,8 @@ def generator( code["name"] = ir.expression.name vals = ", ".join("1" if i else "0" for i in ir.enabled_coefficients) - code["enabled_coefficients"] = f"{{{vals}}}" + # Set to {0} to avoid empty array declaration + code["enabled_coefficients"] = f"{{{vals}}}" if vals != "" else "{{0}}" code["needs_facet_permutations"] = ( "true" if ir.expression.needs_facet_permutations else "false" ) From df8fb7a64ccd7f9ca8b45327108998cc36af4f2b Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 25 Aug 2026 21:50:46 +0200 Subject: [PATCH 08/11] Brackets --- ffcx_backends/cpp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ffcx_backends/cpp.py b/ffcx_backends/cpp.py index 59b18a8..11c7a5f 100644 --- a/ffcx_backends/cpp.py +++ b/ffcx_backends/cpp.py @@ -459,7 +459,7 @@ def generator( vals = ", ".join("1" if i else "0" for i in ir.enabled_coefficients) # Set to {0} to avoid empty array declaration - code["enabled_coefficients"] = f"{{{vals}}}" if vals != "" else "{{0}}" + code["enabled_coefficients"] = f"{{{vals}}}" if vals != "" else "{0}" code["needs_facet_permutations"] = ( "true" if ir.expression.needs_facet_permutations else "false" ) From 8109bca3743bc72efcf20550cabffea47b9cf3cd Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Tue, 25 Aug 2026 22:35:40 +0200 Subject: [PATCH 09/11] Use std::array --- ffcx_backends/cpp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ffcx_backends/cpp.py b/ffcx_backends/cpp.py index 11c7a5f..17853ae 100644 --- a/ffcx_backends/cpp.py +++ b/ffcx_backends/cpp.py @@ -420,7 +420,7 @@ class {factory_name} }} // Data - static constexpr int enabled_coefficients[] = {enabled_coefficients_init}; + static constexpr std::array enabled_coefficients = {enabled_coefficients_init}; static constexpr bool needs_facet_permutations = {needs_facet_permutations}; }}; From fd5ef6c74b95c794f8bb05fe99681ce7bc8d373b Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Wed, 26 Aug 2026 06:39:29 +0200 Subject: [PATCH 10/11] Explicit array template --- ffcx_backends/cpp.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/ffcx_backends/cpp.py b/ffcx_backends/cpp.py index 17853ae..cc837b4 100644 --- a/ffcx_backends/cpp.py +++ b/ffcx_backends/cpp.py @@ -420,7 +420,8 @@ class {factory_name} }} // Data - static constexpr std::array enabled_coefficients = {enabled_coefficients_init}; + static constexpr std::array + enabled_coefficients = {enabled_coefficients_init}; static constexpr bool needs_facet_permutations = {needs_facet_permutations}; }}; @@ -458,8 +459,8 @@ def generator( code["name"] = ir.expression.name vals = ", ".join("1" if i else "0" for i in ir.enabled_coefficients) - # Set to {0} to avoid empty array declaration - code["enabled_coefficients"] = f"{{{vals}}}" if vals != "" else "{0}" + code["num_enabled_coefficients"] = len(ir.enabled_coefficients) + code["enabled_coefficients"] = f"{{{vals}}}" code["needs_facet_permutations"] = ( "true" if ir.expression.needs_facet_permutations else "false" ) @@ -470,6 +471,7 @@ def generator( implementation = integral.factory.format( factory_name=factory_name, enabled_coefficients_init=code["enabled_coefficients"], + num_enabled_coefficients=code["num_enabled_coefficients"], tabulate_tensor=code["tabulate_tensor"], needs_facet_permutations=code["needs_facet_permutations"], scalar_type=options["scalar_type"], From 77f8f77abce1dac1979ee45e8dd4acdf5b16e449 Mon Sep 17 00:00:00 2001 From: Michal Habera Date: Wed, 26 Aug 2026 06:50:37 +0200 Subject: [PATCH 11/11] Mypy --- ffcx_backends/cpp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ffcx_backends/cpp.py b/ffcx_backends/cpp.py index cc837b4..a717ac9 100644 --- a/ffcx_backends/cpp.py +++ b/ffcx_backends/cpp.py @@ -459,7 +459,7 @@ def generator( code["name"] = ir.expression.name vals = ", ".join("1" if i else "0" for i in ir.enabled_coefficients) - code["num_enabled_coefficients"] = len(ir.enabled_coefficients) + code["num_enabled_coefficients"] = str(len(ir.enabled_coefficients)) code["enabled_coefficients"] = f"{{{vals}}}" code["needs_facet_permutations"] = ( "true" if ir.expression.needs_facet_permutations else "false"