Skip to content

Use geometry type independently#4275

Open
michalhabera wants to merge 24 commits into
mainfrom
michal/geom-type
Open

Use geometry type independently#4275
michalhabera wants to merge 24 commits into
mainfrom
michal/geom-type

Conversation

@michalhabera

@michalhabera michalhabera commented Jul 8, 2026

Copy link
Copy Markdown
Contributor

Current use of geometry type different from Function type is half-baked. This PR uses independent type for geometry consistently.

In some places, there is the default behavior to deduce the geometry type as scalar_value_t<T> from the Function type T.

Remember there are two concepts:

  1. dolfinx::scalar T, for Function-like data, could be std::complex<std::float64_t>, and
  2. std::floating_point U for geometry data, cannot be a complex type.

AI assistance: I used Claude Code Opus 4.8 to draft parts of this PR. I reviewed, edited, tested, and take responsibility for the final contribution.

@michalhabera
michalhabera marked this pull request as draft July 8, 2026 14:53
@garth-wells

Copy link
Copy Markdown
Member

It is more subtle than described. The UFCx interface supports a single real scalar type. There is no point mixing, for example, in a straightforward way a double finite element and float geometry because the the float get promoted to double in multiplication, etc. It wasn't half-baked; it was deliberate (but there are of course other ways too).

FEkernel is now more general that the FFCx interface. Is the case where a 'double' FFCx kernel is executed over a float mesh still supported by the PR? Looks like it might not be.

@michalhabera

Copy link
Copy Markdown
Contributor Author

It is more subtle than described. The UFCx interface supports a single real scalar type. There is no point mixing, for example, in a straightforward way a double finite element and float geometry because the the float get promoted to double in multiplication, etc. It wasn't half-baked; it was deliberate (but there are of course other ways too).

FEkernel is now more general that the FFCx interface. Is the case where a 'double' FFCx kernel is executed over a float mesh still supported by the PR? Looks like it might not be.

Right. But why should we limit ourselves to what UFCx/FFCx support at this moment? We've been allowing custom kernels for years.

This is still a draft and I am testing downstream for a combination of scalar_t = float16, geom_t = float64. Think of: I want dolfinx to store a precise float64 mesh, do the assembly loop, but I'll supply my own kernel that does the magic and handles the mixing of geometry of FE fields. And we can improve the UFCx/FFCx do in these cases in near future.

Comment thread python/dolfinx/fem/typemap.py Outdated

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks like nice work towards #3543!

I do not like that this is a global registry: any change in a type needs to be registered not only in the export/wrapper, but also here. Maybe move the lookup table as static member into the python wrapper class? There any newly exported type needs to be made available anyhow. So for Function for example add static Function.cpp_registry, then a general (centralised, but stateless) cpp_type(python_type, dtype) can be implemented that does the lookup (based on the static registry attached to python_type).

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point, I've moved the registry to be local to each class.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is the scalar- and geometry type tuple the right general abstraction? Our exported types currently do not differentiate this, say for Function we have the export Function_float64. Also some types will never depend on such a concept, say for example Vector. I had in mind to just do a dtype based lookup. In case this is not flexible enough (in the future), specialised implementations for these types could be provided.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't fully understand what do you mean. No, this is not a general abstraction, and see my comment #4275 (comment) Some classes will need more, some less types.

With this PR this is possible, but the naming convention is up to you,

dolfinx_wrappers::declare_objects<float, double>(m, "float32_float64");
dolfinx_wrappers::declare_la_objects<std::float16_t>(m, "float16");

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Then the interface (key type) to the registry is not what we want. How about tuple[dtype, ...]? Any export is associated with a variadic length list of types associated to the template types in c++. That we can use as a general interface for type look up.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I've moved registries and respective registry functions to be all local to classes, classmethods. So they have the right signature now, based on how many dtypes they require.

@michalhabera
michalhabera marked this pull request as ready for review July 14, 2026 13:49
@michalhabera

michalhabera commented Jul 14, 2026

Copy link
Copy Markdown
Contributor Author

Updated now, ready for review I think. This is a minimal change and it allows independent use of scalar type and geometry types. I tested e.g. scalar type fp16 and geometry type fp64 downstream.

This is not the whole story, as @garth-wells says. In the next PR I'd make also the kernel output tensor type independent,

std::function<void(V*, const T*, const T*, const U*, const int*, const std::uint8_t*, void*)>

which is the more desired approach, e.g. inputs in low precision and accumulate into higher precision.

@michalhabera
michalhabera requested a review from garth-wells July 14, 2026 16:45
Comment thread python/dolfinx/fem/forms.py Outdated
(np.complex64, np.float32, "complex64"),
(np.complex128, np.float64, "complex128"),
):
register_cpp_type(Form, _scalar, _geom, getattr(_cpp.fem, f"Form_{_name}"))

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just write those out in the class member dict directly? The typing hints of _cpp_object anyhow list them already.

@schnellerhase schnellerhase left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should add a test case for these new templates. Could add a cpp unit test or extend the custom_assembler demo.

@michalhabera

Copy link
Copy Markdown
Contributor Author

We should add a test case for these new templates. Could add a cpp unit test or extend the custom_assembler demo.

I don't plan to do it in this PR. Dolfinx doesn't do anything that would require testing these mixed precision types. And matched precisions are exported and tested. Once this is more explored in derived libraries, we can move some exports incl. their tests to dolfinx.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants