From 80c771b5e9b13cc24db64c79c343492b1bfe9b4e Mon Sep 17 00:00:00 2001 From: Hardik Kothari Date: Tue, 21 Jul 2026 10:54:35 +0200 Subject: [PATCH 1/6] Calibrate symbolic tabulation tolerance from the caller's actual precision --- FIAT/expansions.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 6cd6fa01..c8bc831d 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -9,7 +9,7 @@ import numpy import math -from FIAT import reference_element, jacobi +from FIAT import precision, reference_element, jacobi from FIAT.precision import calibrate_tolerance @@ -379,7 +379,9 @@ def _tabulate(self, n, pts, order=0): if pts.dtype == object: # If binning is undefined, scale by the characteristic function of each subcell - tol = calibrate_tolerance(1E-12, numpy.array(self.ref_el.vertices).dtype) + ref_dtype = numpy.array(self.ref_el.vertices).dtype + dtype = ref_dtype if ref_dtype == numpy.dtype(numpy.float32) else precision.DEFAULT_SCALAR_DTYPE + tol = calibrate_tolerance(1E-12, dtype) Xi = compute_partition_of_unity(self.ref_el, pts, unique=unique, tol=tol) for cell, phi in phis.items(): for alpha in phi: From dcdb65ba7dfdc2b91fe8b73a70cf6882bcf53b41 Mon Sep 17 00:00:00 2001 From: Hardik Kothari Date: Thu, 6 Aug 2026 17:23:15 +0200 Subject: [PATCH 2/6] Propagate dtype through FInAT element construction --- FIAT/expansions.py | 6 ++---- finat/element_factory.py | 25 +++++++++++++++---------- test/finat/test_create_finat_element.py | 23 +++++++++++++++++++++++ 3 files changed, 40 insertions(+), 14 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index c8bc831d..6cd6fa01 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -9,7 +9,7 @@ import numpy import math -from FIAT import precision, reference_element, jacobi +from FIAT import reference_element, jacobi from FIAT.precision import calibrate_tolerance @@ -379,9 +379,7 @@ def _tabulate(self, n, pts, order=0): if pts.dtype == object: # If binning is undefined, scale by the characteristic function of each subcell - ref_dtype = numpy.array(self.ref_el.vertices).dtype - dtype = ref_dtype if ref_dtype == numpy.dtype(numpy.float32) else precision.DEFAULT_SCALAR_DTYPE - tol = calibrate_tolerance(1E-12, dtype) + tol = calibrate_tolerance(1E-12, numpy.array(self.ref_el.vertices).dtype) Xi = compute_partition_of_unity(self.ref_el, pts, unique=unique, tol=tol) for cell, phi in phis.items(): for alpha in phi: diff --git a/finat/element_factory.py b/finat/element_factory.py index 8563aa60..11a7b56a 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -111,13 +111,15 @@ @cache -def as_fiat_cell(cell): +def as_fiat_cell(cell, dtype=None): """Convert a ufl cell to a FIAT cell. - :arg cell: the :class:`ufl.Cell` to convert.""" + :arg cell: the :class:`ufl.Cell` to convert. + :arg dtype: the dtype to build the reference cell's coordinates with. + Defaults to FIAT's own default (float64) if not given.""" if not isinstance(cell, ufl.AbstractCell): raise ValueError("Expecting a UFL Cell") - return ufc_cell(cell) + return ufc_cell(cell, dtype=dtype) @singledispatch @@ -153,7 +155,7 @@ def convert(element, **kwargs): # Base finite elements first @convert.register(finat.ufl.FiniteElement) def convert_finiteelement(element, **kwargs): - cell = as_fiat_cell(element.cell) + cell = as_fiat_cell(element.cell, dtype=kwargs.get("dtype")) if element.family() in {"Quadrature", "Boundary Quadrature"}: degree = element.degree() scheme = element.quadrature_scheme() or "default" @@ -161,7 +163,7 @@ def convert_finiteelement(element, **kwargs): raise ValueError("Quadrature scheme and degree must be specified!") codim = 1 if element.family() == "Boundary Quadrature" else 0 - return finat.make_quadrature_element(cell, degree, scheme, codim), set() + return finat.make_quadrature_element(cell, degree, scheme, codim), {"dtype"} make_finat_element = supported_elements[element.family()] @@ -186,7 +188,7 @@ def convert_finiteelement(element, **kwargs): finat_elem, deps = _create_element(element, **kwargs) return finat.FlattenedDimensions(finat_elem), deps - deps = set() + deps = {"dtype"} finat_kwargs = {} kind = element.variant() if kind is None: @@ -205,7 +207,7 @@ def convert_finiteelement(element, **kwargs): finat_kwargs["variant"] = kind finat_kwargs["shift_axes"] = kwargs["shift_axes"] finat_kwargs["restriction"] = kwargs["restriction"] - deps = {"shift_axes", "restriction"} + deps |= {"shift_axes", "restriction"} else: # Let FIAT handle the general case make_finat_element = finat.Lagrange @@ -227,7 +229,7 @@ def convert_finiteelement(element, **kwargs): finat_kwargs["shift_axes"] = kwargs["shift_axes"] finat_kwargs["restriction"] = kwargs["restriction"] finat_kwargs["continuous"] = False - deps = {"shift_axes", "restriction"} + deps |= {"shift_axes", "restriction"} else: # Let FIAT handle the general case make_finat_element = finat.DiscontinuousLagrange @@ -330,18 +332,21 @@ def convert_restrictedelement(element, **kwargs): _cache = weakref.WeakKeyDictionary() -def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction=None): +def create_element(ufl_element, shape_innermost=True, shift_axes=0, restriction=None, dtype=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) + :arg dtype: the dtype to build the element's reference cell with. + Defaults to FIAT's own default (float64) if not given. """ finat_element, deps = _create_element(ufl_element, shape_innermost=shape_innermost, shift_axes=shift_axes, - restriction=restriction) + restriction=restriction, + dtype=dtype) return finat_element diff --git a/test/finat/test_create_finat_element.py b/test/finat/test_create_finat_element.py index 06113d2b..66624713 100644 --- a/test/finat/test_create_finat_element.py +++ b/test/finat/test_create_finat_element.py @@ -1,3 +1,4 @@ +import numpy import pytest import ufl @@ -166,6 +167,28 @@ def test_cache_hit_vector(ufl_vector_element): assert A is B +def test_dtype_reaches_reference_cell(ufl_element): + """dtype passed to create_element must reach the constructed element's + reference cell, not silently fall back to FIAT's float64 default.""" + default = create_element(ufl_element) + single = create_element(ufl_element, dtype=numpy.float32) + double = create_element(ufl_element, dtype=numpy.float64) + + assert numpy.array(default.cell.vertices).dtype == numpy.float64 + assert numpy.array(single.cell.vertices).dtype == numpy.float32 + assert numpy.array(double.cell.vertices).dtype == numpy.float64 + + +def test_dtype_cache_distinguishes(ufl_element): + """Different dtypes for the same UFL element must not share a cache entry.""" + single_a = create_element(ufl_element, dtype=numpy.float32) + single_b = create_element(ufl_element, dtype=numpy.float32) + double = create_element(ufl_element, dtype=numpy.float64) + + assert single_a is single_b + assert single_a is not double + + if __name__ == "__main__": import os import sys From 6ec3724d98e8930c15dbd75f6706b8794db1ef34 Mon Sep 17 00:00:00 2001 From: Hardik Kothari Date: Fri, 7 Aug 2026 17:03:17 +0200 Subject: [PATCH 3/6] Force double precision for construction-time nullspace/normal computations Rank- and nullspace-determining linear algebra during macro-element construction (compute_normal, spanning_basis, AlfeldSorokinaSpace, hdiv_conforming_coefficients) needs float64 regardless of the working precision, since float32 round-off in vertex coordinates can flip the computed rank. --- FIAT/alfeld_sorokina.py | 48 ++++++++++++++++++++++-------------- FIAT/macro.py | 51 +++++++++++++++++++++++---------------- FIAT/polynomial_set.py | 6 ++++- FIAT/reference_element.py | 5 ++-- 4 files changed, 68 insertions(+), 42 deletions(-) diff --git a/FIAT/alfeld_sorokina.py b/FIAT/alfeld_sorokina.py index e09f26d1..c442e9a4 100644 --- a/FIAT/alfeld_sorokina.py +++ b/FIAT/alfeld_sorokina.py @@ -6,10 +6,13 @@ # # Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024 +import copy + from FIAT import finite_element, dual_set, polynomial_set from FIAT.functional import ComponentPointEvaluation, PointDivergence from FIAT.quadrature_schemes import create_quadrature from FIAT.macro import CkPolynomialSet, AlfeldSplit +from FIAT.reference_element import cast_vertices import numpy @@ -24,25 +27,34 @@ def AlfeldSorokinaSpace(ref_el, degree): num_members = C0.get_num_members() coeffs = C0.get_coeffs() - facet_el = ref_complex.construct_subelement(sd-1) - phi = polynomial_set.ONPolynomialSet(facet_el, 0 if sd == 1 else degree-1) - Q = create_quadrature(facet_el, 2 * phi.degree) - qpts, qwts = Q.get_points(), Q.get_weights() - phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)] - weights = numpy.multiply(phi_at_qpts, qwts) - - rows = [] - for facet in ref_complex.get_interior_facets(sd-1): - n = ref_complex.compute_normal(facet) - jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=1) - div_jump = n[:, None, None] * jumps[1][None, ...] - r = numpy.tensordot(div_jump, weights, axes=(-1, -1)) - rows.append(r.reshape(num_members, -1).T) - - if len(rows) > 0: - dual_mat = numpy.vstack(rows) + interior_facets = ref_complex.get_interior_facets(sd-1) + if len(interior_facets) > 0: + # Redo this in double precision, on a copy of the actual geometry. + ref_el_fp64 = copy.copy(ref_el) + ref_el_fp64.vertices = cast_vertices(ref_el.vertices, float) + ref_el_fp64._split_cache = {} + ref_complex_fp64 = AlfeldSplit(ref_el_fp64) + C0_fp64 = CkPolynomialSet(ref_complex_fp64, degree, order=0, shape=(sd,), variant="bubble") + expansion_set_fp64 = C0_fp64.get_expansion_set() + + facet_el_fp64 = ref_complex_fp64.construct_subelement(sd-1) + phi_fp64 = polynomial_set.ONPolynomialSet(facet_el_fp64, 0 if sd == 1 else degree-1) + Q_fp64 = create_quadrature(facet_el_fp64, 2 * phi_fp64.degree) + qpts_fp64, qwts_fp64 = Q_fp64.get_points(), Q_fp64.get_weights() + phi_at_qpts_fp64 = phi_fp64.tabulate(qpts_fp64)[(0,) * (sd-1)] + weights_fp64 = numpy.multiply(phi_at_qpts_fp64, qwts_fp64) + + rows_fp64 = [] + for facet in ref_complex_fp64.get_interior_facets(sd-1): + n_fp64 = ref_complex_fp64.compute_normal(facet) + jumps_fp64 = expansion_set_fp64.tabulate_normal_jumps(degree, qpts_fp64, facet, order=1) + div_jump_fp64 = n_fp64[:, None, None] * jumps_fp64[1][None, ...] + r_fp64 = numpy.tensordot(div_jump_fp64, weights_fp64, axes=(-1, -1)) + rows_fp64.append(r_fp64.reshape(num_members, -1).T) + + dual_mat = numpy.vstack(rows_fp64) nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) - coeffs = numpy.tensordot(nsp, coeffs, axes=(-1, 0)) + coeffs = numpy.tensordot(nsp.astype(coeffs.dtype), coeffs, axes=(-1, 0)) return polynomial_set.PolynomialSet(ref_complex, degree, degree, expansion_set, coeffs) diff --git a/FIAT/macro.py b/FIAT/macro.py index 2f2795b1..637accd0 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -1,3 +1,4 @@ +import copy from itertools import chain, combinations import numpy @@ -528,29 +529,37 @@ def hdiv_conforming_coefficients(U, order=0): k = 1 if expansion_set.continuity == "C0" else 0 sd = ref_el.get_spatial_dimension() - facet_el = ref_el.construct_subelement(sd-1) - phi_deg = 0 if sd == 1 else degree - k - phi = polynomial_set.ONPolynomialSet(facet_el, phi_deg, shape=shape[1:]) - Q = create_quadrature(facet_el, 2 * phi_deg) - qpts, qwts = Q.get_points(), Q.get_weights() - phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)] - weights = numpy.multiply(phi_at_qpts, qwts) - ax = tuple(range(1, weights.ndim)) - - rows = [] - for facet in ref_el.get_interior_facets(sd-1): - normal = ref_el.compute_scaled_normal(facet) - ncoeffs = numpy.tensordot(coeffs, normal, axes=(len(shape), 0)) - jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=order) - for r in range(k, order+1): - njump = numpy.dot(ncoeffs, jumps[r]) - rows.append(numpy.tensordot(weights, njump, axes=(ax, ax))) - - if len(rows) > 0: - dual_mat = numpy.vstack(rows) + + interior_facets = ref_el.get_interior_facets(sd-1) + if len(interior_facets) > 0: + # Redo this in double precision, on a copy of the actual geometry. + parent_fp64 = copy.copy(ref_el.get_parent()) + parent_fp64.vertices = reference_element.cast_vertices(parent_fp64.vertices, float) + parent_fp64._split_cache = {} + ref_complex_fp64 = type(ref_el)(parent_fp64) + expansion_set_fp64 = expansions.ExpansionSet(ref_complex_fp64, scale=expansion_set.scale, variant=expansion_set.variant) + facet_el_fp64 = ref_complex_fp64.construct_subelement(sd-1) + phi_fp64 = polynomial_set.ONPolynomialSet(facet_el_fp64, phi_deg, shape=shape[1:]) + Q_fp64 = create_quadrature(facet_el_fp64, 2 * phi_deg) + qpts_fp64, qwts_fp64 = Q_fp64.get_points(), Q_fp64.get_weights() + phi_at_qpts_fp64 = phi_fp64.tabulate(qpts_fp64)[(0,) * (sd-1)] + weights_fp64 = numpy.multiply(phi_at_qpts_fp64, qwts_fp64) + ax = tuple(range(1, weights_fp64.ndim)) + coeffs_fp64 = coeffs.astype(numpy.float64) + + rows_fp64 = [] + for facet in ref_complex_fp64.get_interior_facets(sd-1): + normal_fp64 = ref_complex_fp64.compute_scaled_normal(facet) + ncoeffs_fp64 = numpy.tensordot(coeffs_fp64, normal_fp64, axes=(len(shape), 0)) + jumps_fp64 = expansion_set_fp64.tabulate_normal_jumps(degree, qpts_fp64, facet, order=order) + for r in range(k, order+1): + njump_fp64 = numpy.dot(ncoeffs_fp64, jumps_fp64[r]) + rows_fp64.append(numpy.tensordot(weights_fp64, njump_fp64, axes=(ax, ax))) + + dual_mat = numpy.vstack(rows_fp64) nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) - coeffs = numpy.tensordot(nsp, coeffs, axes=(1, 0)) + coeffs = numpy.tensordot(nsp.astype(coeffs.dtype), coeffs, axes=(1, 0)) return coeffs diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index b3aee67f..c42b298c 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -16,6 +16,8 @@ # an entire set of polynomials) import numpy + +from FIAT.precision import calibrate_tolerance from itertools import chain from FIAT import expansions @@ -157,9 +159,11 @@ def form_matrix_product(mats, alpha): return result -def spanning_basis(A, nullspace=False, rtol=1e-10): +def spanning_basis(A, nullspace=False, rtol=None): """Construct a basis that spans the rows of A via SVD. """ + if rtol is None: + rtol = calibrate_tolerance(1e-10, A.dtype) Aflat = A.reshape(A.shape[0], -1) u, sig, vt = numpy.linalg.svd(Aflat, full_matrices=True) atol = rtol * (sig[0] + 1) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 3427e93c..962c8b28 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -409,7 +409,8 @@ def compute_normal(self, facet_i, cell=None): if cell is None: cell = next(k for k, facets in enumerate(self.connectivity[(sd, sd-1)]) if facet_i in facets) - verts = numpy.asarray(self.get_vertices_of_subcomplex(t[sd][cell])) + # Always compute in double precision for a robust SVD below. + verts = numpy.asarray(self.get_vertices_of_subcomplex(t[sd][cell]), dtype=float) # Interval case if self.get_shape() == LINE: v_i = t[1][cell].index(t[0][facet_i][0]) @@ -429,7 +430,7 @@ def compute_normal(self, facet_i, cell=None): self.get_vertices_of_subcomplex(t[sd-1][facet_i]) # now I find everything normal to the facet. - vcf = numpy.asarray(vert_coords_of_facet) + vcf = numpy.asarray(vert_coords_of_facet, dtype=float) facet_span = vcf[1:, :] - vcf[:1, :] (_, sf, vft) = numpy.linalg.svd(facet_span) From 14f55f79acdfdc37497e98b40a4f82cd0213ed30 Mon Sep 17 00:00:00 2001 From: Hardik Kothari Date: Thu, 3 Sep 2026 16:41:38 +0200 Subject: [PATCH 4/6] Fix fp32 macro-element rank construction --- FIAT/macro.py | 112 +++++++++++++++++++++++--- test/FIAT/unit/test_hct.py | 15 ++-- test/FIAT/unit/test_macro.py | 13 ++- test/FIAT/unit/test_powell_sabin.py | 10 ++- test/FIAT/unit/test_stokes_complex.py | 17 +++- 5 files changed, 139 insertions(+), 28 deletions(-) diff --git a/FIAT/macro.py b/FIAT/macro.py index 08298aef..e78e15d5 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -88,7 +88,10 @@ class SplitSimplicialComplex(SimplicialComplex): :arg vertices: The vertices of the simplicial complex. :arg topology: The topology of the simplicial complex. """ - def __init__(self, parent, vertices, topology): + def __init__( + self, parent: SimplicialComplex, vertices: tuple, topology: dict, + split_parent: SimplicialComplex | None = None) -> None: + self._split_parent = parent if split_parent is None else split_parent self._parent_complex = parent while parent.get_parent(): parent = parent.get_parent() @@ -199,6 +202,22 @@ def get_parent(self): def get_parent_complex(self): return self._parent_complex + def reconstruct( + self, ref_el: SimplicialComplex) -> "SplitSimplicialComplex": + """Reconstruct this split on another reference complex. + + Parameters + ---------- + ref_el : SimplicialComplex + The reference complex to split. + + Returns + ------- + SplitSimplicialComplex + The reconstructed split. + """ + return type(self)(ref_el) + class IsoSplit(SplitSimplicialComplex): """Splits simplex into the simplicial complex obtained by @@ -237,6 +256,21 @@ def __init__(self, ref_el, degree=2, variant=None): new_topology = make_topology(sd, len(new_verts), edges) super().__init__(ref_el, tuple(new_verts), new_topology) + def reconstruct(self, ref_el: SimplicialComplex) -> "IsoSplit": + """Reconstruct this split on another reference complex. + + Parameters + ---------- + ref_el : SimplicialComplex + The reference complex to split. + + Returns + ------- + IsoSplit + The reconstructed split. + """ + return type(self)(ref_el, degree=self.degree, variant=self.variant) + def construct_subcomplex(self, dimension): """Constructs the reference subcomplex of the parent complex specified by subcomplex dimension. @@ -288,7 +322,25 @@ def __init__(self, ref_el, dimension=1): new_topology[sd] = dict(enumerate(simplices)) parent = ref_el if dimension == sd else PowellSabinSplit(ref_el, dimension=dimension+1) - super().__init__(parent, tuple(new_verts), new_topology) + super().__init__(parent, tuple(new_verts), new_topology, + split_parent=ref_el) + + def reconstruct(self, ref_el: SimplicialComplex) -> "PowellSabinSplit": + """Reconstruct this split on another reference complex. + + Parameters + ---------- + ref_el : SimplicialComplex + The reference complex to split. + + Returns + ------- + PowellSabinSplit + The reconstructed split. + """ + if type(self) is PowellSabinSplit: + return type(self)(ref_el, dimension=self.split_dimension) + return super().reconstruct(ref_el) def construct_subcomplex(self, dimension): """Constructs the reference subcomplex of the parent complex @@ -345,15 +397,16 @@ def __init__(self, ref_el): assert ref_el.get_shape() == TRIANGLE verts = ref_el.get_vertices() new_verts = list(verts) + dtype = numpy.asarray(verts).dtype new_verts.extend( map(tuple, bary_to_xy(verts, - [(1/3, 1/3, 1/3), - (1/2, 1/2, 0), - (1/2, 0, 1/2), - (0, 1/2, 1/2), - (1/2, 1/4, 1/4), - (1/4, 1/2, 1/4), - (1/4, 1/4, 1/2)]))) + numpy.asarray([(1/3, 1/3, 1/3), + (1/2, 1/2, 0), + (1/2, 0, 1/2), + (0, 1/2, 1/2), + (1/2, 1/4, 1/4), + (1/4, 1/2, 1/4), + (1/4, 1/4, 1/2)], dtype=dtype)))) edges = [(0, 4), (0, 7), (0, 5), (1, 4), (1, 8), (1, 6), @@ -363,7 +416,8 @@ def __init__(self, ref_el): parent = PowellSabinSplit(ref_el) new_topology = make_topology(2, len(new_verts), edges) - super().__init__(parent, tuple(new_verts), new_topology) + super().__init__(parent, tuple(new_verts), new_topology, + split_parent=ref_el) def construct_subcomplex(self, dimension): """Constructs the reference subcomplex of the parent cell subentity @@ -379,6 +433,31 @@ def construct_subcomplex(self, dimension): raise ValueError("Illegal dimension") +def _reconstruct_split_complex_fp64( + ref_el: SimplicialComplex) -> SimplicialComplex: + """Reconstruct a split complex from float64 root geometry. + + Parameters + ---------- + ref_el : SimplicialComplex + The reference complex to reconstruct. + + Returns + ------- + SimplicialComplex + The reconstructed reference complex. + """ + if not isinstance(ref_el, SplitSimplicialComplex): + ref_el_fp64 = copy.copy(ref_el) + ref_el_fp64.vertices = reference_element.cast_vertices( + ref_el.vertices, numpy.float64) + ref_el_fp64._split_cache = {} + return ref_el_fp64 + + parent_fp64 = _reconstruct_split_complex_fp64(ref_el._split_parent) + return ref_el.reconstruct(parent_fp64) + + class MacroQuadratureRule(QuadratureRule): """Composite quadrature rule on parent facets that respects the splitting. @@ -451,6 +530,18 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs): if not isinstance(order, (int, dict)): raise TypeError(f"'order' must be either an int or dict, not {type(order).__name__}") + expansion_set = expansions.ExpansionSet(ref_el, **kwargs) + dtype = numpy.asarray(ref_el.get_vertices()).dtype + if dtype == numpy.dtype(numpy.float32): + ref_el_fp64 = _reconstruct_split_complex_fp64(ref_el) + poly_set_fp64 = CkPolynomialSet( + ref_el_fp64, degree, order=copy.deepcopy(order), + vorder=vorder, shape=shape, **kwargs) + coeffs = poly_set_fp64.get_coeffs() + super().__init__( + ref_el, degree, degree, expansion_set, coeffs) + return + sd = ref_el.get_spatial_dimension() if isinstance(order, int): order = {sd-1: dict.fromkeys(ref_el.get_interior_facets(sd-1), order)} @@ -462,7 +553,6 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs): if not all(k in {0, sd-1} for k in order): raise NotImplementedError("Only face or vertex constraints have been implemented.") - expansion_set = expansions.ExpansionSet(ref_el, **kwargs) k = 1 if expansion_set.continuity == "C0" else 0 # Impose C^forder continuity across interior facets diff --git a/test/FIAT/unit/test_hct.py b/test/FIAT/unit/test_hct.py index cb2d56fb..73109f91 100644 --- a/test/FIAT/unit/test_hct.py +++ b/test/FIAT/unit/test_hct.py @@ -7,10 +7,11 @@ from FIAT.macro import CkPolynomialSet -@pytest.fixture -def cell(): - K = ufc_simplex(2) - K.vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84)) +@pytest.fixture(params=(numpy.float64, numpy.float32), ids=("float64", "float32")) +def cell(request): + K = ufc_simplex(2, dtype=request.param) + K.vertices = tuple(map(tuple, numpy.asarray( + ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84)), dtype=request.param))) return K @@ -76,7 +77,9 @@ def test_full_polynomials(cell, reduced): assert span_greater_equal(tab, C1_tab) +@pytest.mark.parametrize("dtype", (numpy.float64, numpy.float32), + ids=("float64", "float32")) @pytest.mark.parametrize("degree, space_dimension", [(13, 127), (14, 144)]) -def test_hct_high_order_degree(degree: int, space_dimension: int) -> None: - fe = HCT(ufc_simplex(2), degree) +def test_hct_high_order_degree(dtype, degree: int, space_dimension: int) -> None: + fe = HCT(ufc_simplex(2, dtype=dtype), degree) assert fe.space_dimension() == space_dimension diff --git a/test/FIAT/unit/test_macro.py b/test/FIAT/unit/test_macro.py index 9ad45a95..db2c6450 100644 --- a/test/FIAT/unit/test_macro.py +++ b/test/FIAT/unit/test_macro.py @@ -375,11 +375,15 @@ def test_macro_expansion(cell, split, variant, degree): @pytest.mark.parametrize("order", (0, 1)) @pytest.mark.parametrize("variant", (None, "bubble")) @pytest.mark.parametrize("degree", (1, 4)) -def test_Ck_basis(cell, order, degree, variant): +@pytest.mark.parametrize("split", (AlfeldSplit, IsoSplit)) +@pytest.mark.parametrize("dtype", (numpy.float64, numpy.float32), + ids=("float64", "float32")) +def test_Ck_basis(cell, order, degree, variant, split, dtype): # Test that we can correctly tabulate on points on facets. # This breaks if we were binning points into more than one cell without a partition of unity. # It suffices to tabulate on the vertices of the simplicial complex. - A = AlfeldSplit(cell) + cell = ufc_simplex(cell.get_spatial_dimension(), dtype=dtype) + A = split(cell) Ck = CkPolynomialSet(A, degree, order=order, variant=variant) U = Ck.get_expansion_set() cell_node_map = U.get_cell_node_map(degree) @@ -397,11 +401,12 @@ def test_Ck_basis(cell, order, degree, variant): assert numpy.allclose(local_phis, phis[:, ipts]) -def test_C2_double_alfeld(): +@pytest.mark.parametrize("dtype", (numpy.float64, numpy.float32)) +def test_C2_double_alfeld(dtype): from FIAT.c2_elements import AlfeldC2Space # Construct the quintic C2 spline on the double Alfeld split # See Section 7.5 of Lai & Schumacher - K = ufc_simplex(2) + K = ufc_simplex(2, dtype=dtype) degree = 5 P = AlfeldC2Space(K, degree) assert P.get_num_members() == 27 diff --git a/test/FIAT/unit/test_powell_sabin.py b/test/FIAT/unit/test_powell_sabin.py index f97c7104..31f19ed7 100644 --- a/test/FIAT/unit/test_powell_sabin.py +++ b/test/FIAT/unit/test_powell_sabin.py @@ -6,15 +6,19 @@ from FIAT.reference_element import make_lattice, ufc_simplex -@pytest.fixture -def cell(): - return ufc_simplex(2) +@pytest.fixture(params=(numpy.float64, numpy.float32), + ids=("float64", "float32")) +def cell(request): + return ufc_simplex(2, dtype=request.param) @pytest.mark.parametrize("el", (PS6, PS12)) def test_powell_sabin_constant(cell, el): # Test that bfs associated with point evaluation sum up to 1 fe = el(cell) + assert (numpy.asarray(fe.get_reference_complex().vertices).dtype + == numpy.asarray(cell.vertices).dtype) + assert fe.get_coeffs().dtype == numpy.float64 pts = make_lattice(cell.get_vertices(), 3) tab = fe.tabulate(2, pts) diff --git a/test/FIAT/unit/test_stokes_complex.py b/test/FIAT/unit/test_stokes_complex.py index a7725ed7..0eb4b46c 100644 --- a/test/FIAT/unit/test_stokes_complex.py +++ b/test/FIAT/unit/test_stokes_complex.py @@ -143,8 +143,10 @@ def check_stokes_complex(spaces, degree): @pytest.mark.parametrize("reduced", (False, True), ids=("full", "reduced")) @pytest.mark.parametrize("sobolev", ("H1", "H1div")) -@pytest.mark.parametrize("cell", (T,)) -def test_hct_stokes_complex(cell, sobolev, reduced): +@pytest.mark.parametrize("dtype", (numpy.float64, numpy.float32), + ids=("float64", "float32")) +def test_hct_stokes_complex(dtype, sobolev, reduced): + cell = ufc_simplex(2, dtype=dtype) if sobolev == "H1": if reduced: spaces = [rHCT(cell), rAQ(cell), DG(cell, 0)] @@ -246,11 +248,18 @@ def test_gn_trace(sd): @pytest.mark.parametrize("cell", (T, S)) -@pytest.mark.parametrize("family", ("AQ", "CH", "GN", "GN2")) -def test_minimal_stokes_space(cell, family): +@pytest.mark.parametrize("family,dtype", ( + ("AQ", numpy.float64), + ("CH", numpy.float64), + ("CH", numpy.float32), + ("GN", numpy.float64), + ("GN2", numpy.float64), +)) +def test_minimal_stokes_space(cell, family, dtype): # Test that the C0 Stokes space is spanned by a C0 basis # Also test that its divergence is constant sd = cell.get_spatial_dimension() + cell = ufc_simplex(sd, dtype=dtype) if family == "GN": degree = 1 space = GuzmanNeilanSpace From 7e2b3326c08226fec3b050b1d82f2bd3325e97db Mon Sep 17 00:00:00 2001 From: Hardik Kothari Date: Fri, 4 Sep 2026 12:20:44 +0200 Subject: [PATCH 5/6] Keep FIAT construction in double precision --- FIAT/alfeld_sorokina.py | 48 +++---- FIAT/expansions.py | 2 +- FIAT/macro.py | 164 +++++------------------- FIAT/polynomial_set.py | 5 +- FIAT/reference_element.py | 23 ++-- test/FIAT/unit/test_hct.py | 3 +- test/FIAT/unit/test_powell_sabin.py | 4 +- test/finat/test_create_finat_element.py | 8 +- 8 files changed, 73 insertions(+), 184 deletions(-) diff --git a/FIAT/alfeld_sorokina.py b/FIAT/alfeld_sorokina.py index c442e9a4..e09f26d1 100644 --- a/FIAT/alfeld_sorokina.py +++ b/FIAT/alfeld_sorokina.py @@ -6,13 +6,10 @@ # # Written by Pablo D. Brubeck (brubeck@protonmail.com), 2024 -import copy - from FIAT import finite_element, dual_set, polynomial_set from FIAT.functional import ComponentPointEvaluation, PointDivergence from FIAT.quadrature_schemes import create_quadrature from FIAT.macro import CkPolynomialSet, AlfeldSplit -from FIAT.reference_element import cast_vertices import numpy @@ -27,34 +24,25 @@ def AlfeldSorokinaSpace(ref_el, degree): num_members = C0.get_num_members() coeffs = C0.get_coeffs() - interior_facets = ref_complex.get_interior_facets(sd-1) - if len(interior_facets) > 0: - # Redo this in double precision, on a copy of the actual geometry. - ref_el_fp64 = copy.copy(ref_el) - ref_el_fp64.vertices = cast_vertices(ref_el.vertices, float) - ref_el_fp64._split_cache = {} - ref_complex_fp64 = AlfeldSplit(ref_el_fp64) - C0_fp64 = CkPolynomialSet(ref_complex_fp64, degree, order=0, shape=(sd,), variant="bubble") - expansion_set_fp64 = C0_fp64.get_expansion_set() - - facet_el_fp64 = ref_complex_fp64.construct_subelement(sd-1) - phi_fp64 = polynomial_set.ONPolynomialSet(facet_el_fp64, 0 if sd == 1 else degree-1) - Q_fp64 = create_quadrature(facet_el_fp64, 2 * phi_fp64.degree) - qpts_fp64, qwts_fp64 = Q_fp64.get_points(), Q_fp64.get_weights() - phi_at_qpts_fp64 = phi_fp64.tabulate(qpts_fp64)[(0,) * (sd-1)] - weights_fp64 = numpy.multiply(phi_at_qpts_fp64, qwts_fp64) - - rows_fp64 = [] - for facet in ref_complex_fp64.get_interior_facets(sd-1): - n_fp64 = ref_complex_fp64.compute_normal(facet) - jumps_fp64 = expansion_set_fp64.tabulate_normal_jumps(degree, qpts_fp64, facet, order=1) - div_jump_fp64 = n_fp64[:, None, None] * jumps_fp64[1][None, ...] - r_fp64 = numpy.tensordot(div_jump_fp64, weights_fp64, axes=(-1, -1)) - rows_fp64.append(r_fp64.reshape(num_members, -1).T) - - dual_mat = numpy.vstack(rows_fp64) + facet_el = ref_complex.construct_subelement(sd-1) + phi = polynomial_set.ONPolynomialSet(facet_el, 0 if sd == 1 else degree-1) + Q = create_quadrature(facet_el, 2 * phi.degree) + qpts, qwts = Q.get_points(), Q.get_weights() + phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)] + weights = numpy.multiply(phi_at_qpts, qwts) + + rows = [] + for facet in ref_complex.get_interior_facets(sd-1): + n = ref_complex.compute_normal(facet) + jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=1) + div_jump = n[:, None, None] * jumps[1][None, ...] + r = numpy.tensordot(div_jump, weights, axes=(-1, -1)) + rows.append(r.reshape(num_members, -1).T) + + if len(rows) > 0: + dual_mat = numpy.vstack(rows) nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) - coeffs = numpy.tensordot(nsp.astype(coeffs.dtype), coeffs, axes=(-1, 0)) + coeffs = numpy.tensordot(nsp, coeffs, axes=(-1, 0)) return polynomial_set.PolynomialSet(ref_complex, degree, degree, expansion_set, coeffs) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 33f3ed00..52f224bb 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -459,7 +459,7 @@ def _tabulate(self, n, pts, order=0): if pts.dtype == object: # If binning is undefined, scale by the characteristic function of each subcell - tol = calibrate_tolerance(1E-12, numpy.array(self.ref_el.vertices).dtype) + tol = calibrate_tolerance(1E-12, self.ref_el.target_dtype) Xi = compute_partition_of_unity(self.ref_el, pts, unique=unique, tol=tol) for cell, phi in phis.items(): for alpha in phi: diff --git a/FIAT/macro.py b/FIAT/macro.py index e78e15d5..36aad384 100644 --- a/FIAT/macro.py +++ b/FIAT/macro.py @@ -1,4 +1,3 @@ -import copy from itertools import chain, combinations import numpy @@ -88,10 +87,7 @@ class SplitSimplicialComplex(SimplicialComplex): :arg vertices: The vertices of the simplicial complex. :arg topology: The topology of the simplicial complex. """ - def __init__( - self, parent: SimplicialComplex, vertices: tuple, topology: dict, - split_parent: SimplicialComplex | None = None) -> None: - self._split_parent = parent if split_parent is None else split_parent + def __init__(self, parent, vertices, topology): self._parent_complex = parent while parent.get_parent(): parent = parent.get_parent() @@ -155,6 +151,7 @@ def __init__( self._interior_facets = interior_facets super().__init__(parent.shape, vertices, topology) + self.target_dtype = parent.target_dtype def get_child_to_parent(self): """Maps split complex facet tuple to its parent entity tuple.""" @@ -202,22 +199,6 @@ def get_parent(self): def get_parent_complex(self): return self._parent_complex - def reconstruct( - self, ref_el: SimplicialComplex) -> "SplitSimplicialComplex": - """Reconstruct this split on another reference complex. - - Parameters - ---------- - ref_el : SimplicialComplex - The reference complex to split. - - Returns - ------- - SplitSimplicialComplex - The reconstructed split. - """ - return type(self)(ref_el) - class IsoSplit(SplitSimplicialComplex): """Splits simplex into the simplicial complex obtained by @@ -256,21 +237,6 @@ def __init__(self, ref_el, degree=2, variant=None): new_topology = make_topology(sd, len(new_verts), edges) super().__init__(ref_el, tuple(new_verts), new_topology) - def reconstruct(self, ref_el: SimplicialComplex) -> "IsoSplit": - """Reconstruct this split on another reference complex. - - Parameters - ---------- - ref_el : SimplicialComplex - The reference complex to split. - - Returns - ------- - IsoSplit - The reconstructed split. - """ - return type(self)(ref_el, degree=self.degree, variant=self.variant) - def construct_subcomplex(self, dimension): """Constructs the reference subcomplex of the parent complex specified by subcomplex dimension. @@ -322,25 +288,7 @@ def __init__(self, ref_el, dimension=1): new_topology[sd] = dict(enumerate(simplices)) parent = ref_el if dimension == sd else PowellSabinSplit(ref_el, dimension=dimension+1) - super().__init__(parent, tuple(new_verts), new_topology, - split_parent=ref_el) - - def reconstruct(self, ref_el: SimplicialComplex) -> "PowellSabinSplit": - """Reconstruct this split on another reference complex. - - Parameters - ---------- - ref_el : SimplicialComplex - The reference complex to split. - - Returns - ------- - PowellSabinSplit - The reconstructed split. - """ - if type(self) is PowellSabinSplit: - return type(self)(ref_el, dimension=self.split_dimension) - return super().reconstruct(ref_el) + super().__init__(parent, tuple(new_verts), new_topology) def construct_subcomplex(self, dimension): """Constructs the reference subcomplex of the parent complex @@ -397,16 +345,15 @@ def __init__(self, ref_el): assert ref_el.get_shape() == TRIANGLE verts = ref_el.get_vertices() new_verts = list(verts) - dtype = numpy.asarray(verts).dtype new_verts.extend( map(tuple, bary_to_xy(verts, - numpy.asarray([(1/3, 1/3, 1/3), - (1/2, 1/2, 0), - (1/2, 0, 1/2), - (0, 1/2, 1/2), - (1/2, 1/4, 1/4), - (1/4, 1/2, 1/4), - (1/4, 1/4, 1/2)], dtype=dtype)))) + [(1/3, 1/3, 1/3), + (1/2, 1/2, 0), + (1/2, 0, 1/2), + (0, 1/2, 1/2), + (1/2, 1/4, 1/4), + (1/4, 1/2, 1/4), + (1/4, 1/4, 1/2)]))) edges = [(0, 4), (0, 7), (0, 5), (1, 4), (1, 8), (1, 6), @@ -416,8 +363,7 @@ def __init__(self, ref_el): parent = PowellSabinSplit(ref_el) new_topology = make_topology(2, len(new_verts), edges) - super().__init__(parent, tuple(new_verts), new_topology, - split_parent=ref_el) + super().__init__(parent, tuple(new_verts), new_topology) def construct_subcomplex(self, dimension): """Constructs the reference subcomplex of the parent cell subentity @@ -433,31 +379,6 @@ def construct_subcomplex(self, dimension): raise ValueError("Illegal dimension") -def _reconstruct_split_complex_fp64( - ref_el: SimplicialComplex) -> SimplicialComplex: - """Reconstruct a split complex from float64 root geometry. - - Parameters - ---------- - ref_el : SimplicialComplex - The reference complex to reconstruct. - - Returns - ------- - SimplicialComplex - The reconstructed reference complex. - """ - if not isinstance(ref_el, SplitSimplicialComplex): - ref_el_fp64 = copy.copy(ref_el) - ref_el_fp64.vertices = reference_element.cast_vertices( - ref_el.vertices, numpy.float64) - ref_el_fp64._split_cache = {} - return ref_el_fp64 - - parent_fp64 = _reconstruct_split_complex_fp64(ref_el._split_parent) - return ref_el.reconstruct(parent_fp64) - - class MacroQuadratureRule(QuadratureRule): """Composite quadrature rule on parent facets that respects the splitting. @@ -530,18 +451,6 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs): if not isinstance(order, (int, dict)): raise TypeError(f"'order' must be either an int or dict, not {type(order).__name__}") - expansion_set = expansions.ExpansionSet(ref_el, **kwargs) - dtype = numpy.asarray(ref_el.get_vertices()).dtype - if dtype == numpy.dtype(numpy.float32): - ref_el_fp64 = _reconstruct_split_complex_fp64(ref_el) - poly_set_fp64 = CkPolynomialSet( - ref_el_fp64, degree, order=copy.deepcopy(order), - vorder=vorder, shape=shape, **kwargs) - coeffs = poly_set_fp64.get_coeffs() - super().__init__( - ref_el, degree, degree, expansion_set, coeffs) - return - sd = ref_el.get_spatial_dimension() if isinstance(order, int): order = {sd-1: dict.fromkeys(ref_el.get_interior_facets(sd-1), order)} @@ -553,6 +462,7 @@ def __init__(self, ref_el, degree, order=1, vorder=None, shape=(), **kwargs): if not all(k in {0, sd-1} for k in order): raise NotImplementedError("Only face or vertex constraints have been implemented.") + expansion_set = expansions.ExpansionSet(ref_el, **kwargs) k = 1 if expansion_set.continuity == "C0" else 0 # Impose C^forder continuity across interior facets @@ -630,37 +540,29 @@ def hdiv_conforming_coefficients(U, order=0): k = 1 if expansion_set.continuity == "C0" else 0 sd = ref_el.get_spatial_dimension() - phi_deg = 0 if sd == 1 else degree - k + facet_el = ref_el.construct_subelement(sd-1) - interior_facets = ref_el.get_interior_facets(sd-1) - if len(interior_facets) > 0: - # Redo this in double precision, on a copy of the actual geometry. - parent_fp64 = copy.copy(ref_el.get_parent()) - parent_fp64.vertices = reference_element.cast_vertices(parent_fp64.vertices, float) - parent_fp64._split_cache = {} - ref_complex_fp64 = type(ref_el)(parent_fp64) - expansion_set_fp64 = expansions.ExpansionSet(ref_complex_fp64, scale=expansion_set.scale, variant=expansion_set.variant) - facet_el_fp64 = ref_complex_fp64.construct_subelement(sd-1) - phi_fp64 = polynomial_set.ONPolynomialSet(facet_el_fp64, phi_deg, shape=shape[1:]) - Q_fp64 = create_quadrature(facet_el_fp64, 2 * phi_deg) - qpts_fp64, qwts_fp64 = Q_fp64.get_points(), Q_fp64.get_weights() - phi_at_qpts_fp64 = phi_fp64.tabulate(qpts_fp64)[(0,) * (sd-1)] - weights_fp64 = numpy.multiply(phi_at_qpts_fp64, qwts_fp64) - ax = tuple(range(1, weights_fp64.ndim)) - coeffs_fp64 = coeffs.astype(numpy.float64) - - rows_fp64 = [] - for facet in ref_complex_fp64.get_interior_facets(sd-1): - normal_fp64 = ref_complex_fp64.compute_scaled_normal(facet) - ncoeffs_fp64 = numpy.tensordot(coeffs_fp64, normal_fp64, axes=(len(shape), 0)) - jumps_fp64 = expansion_set_fp64.tabulate_normal_jumps(degree, qpts_fp64, facet, order=order) - for r in range(k, order+1): - njump_fp64 = numpy.dot(ncoeffs_fp64, jumps_fp64[r]) - rows_fp64.append(numpy.tensordot(weights_fp64, njump_fp64, axes=(ax, ax))) - - dual_mat = numpy.vstack(rows_fp64) + phi_deg = 0 if sd == 1 else degree - k + phi = polynomial_set.ONPolynomialSet(facet_el, phi_deg, shape=shape[1:]) + Q = create_quadrature(facet_el, 2 * phi_deg) + qpts, qwts = Q.get_points(), Q.get_weights() + phi_at_qpts = phi.tabulate(qpts)[(0,) * (sd-1)] + weights = numpy.multiply(phi_at_qpts, qwts) + ax = tuple(range(1, weights.ndim)) + + rows = [] + for facet in ref_el.get_interior_facets(sd-1): + normal = ref_el.compute_scaled_normal(facet) + ncoeffs = numpy.tensordot(coeffs, normal, axes=(len(shape), 0)) + jumps = expansion_set.tabulate_normal_jumps(degree, qpts, facet, order=order) + for r in range(k, order+1): + njump = numpy.dot(ncoeffs, jumps[r]) + rows.append(numpy.tensordot(weights, njump, axes=(ax, ax))) + + if len(rows) > 0: + dual_mat = numpy.vstack(rows) nsp = polynomial_set.spanning_basis(dual_mat, nullspace=True) - coeffs = numpy.tensordot(nsp.astype(coeffs.dtype), coeffs, axes=(1, 0)) + coeffs = numpy.tensordot(nsp, coeffs, axes=(1, 0)) return coeffs diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index c42b298c..4b518832 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -17,7 +17,6 @@ import numpy -from FIAT.precision import calibrate_tolerance from itertools import chain from FIAT import expansions @@ -159,11 +158,9 @@ def form_matrix_product(mats, alpha): return result -def spanning_basis(A, nullspace=False, rtol=None): +def spanning_basis(A, nullspace=False, rtol=1e-10): """Construct a basis that spans the rows of A via SVD. """ - if rtol is None: - rtol = calibrate_tolerance(1e-10, A.dtype) Aflat = A.reshape(A.shape[0], -1) u, sig, vt = numpy.linalg.svd(Aflat, full_matrices=True) atol = rtol * (sig[0] + 1) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 962c8b28..829debb1 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -143,6 +143,7 @@ def __init__(self, shape, vertices, topology): comprising the facet.""" self.shape = shape self.vertices = vertices + self.target_dtype = numpy.asarray(vertices).dtype self.topology = topology # Given the topology, work out for each entity in the cell, @@ -1659,9 +1660,7 @@ def ufc_hypercube(spatial_dim, dtype=None): """Factory function that maps spatial dimension to an instance of the UFC reference hypercube of that dimension. - :arg dtype: optional numpy dtype to cast the vertex coordinates to. - Defaults to the plain Python `float` used by the hardcoded - vertex coordinates. + :arg dtype: optional working dtype for tabulation. """ if spatial_dim == 0: cell = Point() @@ -1674,7 +1673,7 @@ def ufc_hypercube(spatial_dim, dtype=None): else: raise RuntimeError(f"Can't create UFC hypercube of dimension {spatial_dim}.") if dtype is not None: - cell.vertices = cast_vertices(cell.vertices, dtype) + cell.target_dtype = numpy.dtype(dtype) return cell @@ -1697,9 +1696,7 @@ def ufc_simplex(spatial_dim, dtype=None): """Factory function that maps spatial dimension to an instance of the UFC reference simplex of that dimension. - :arg dtype: optional numpy dtype to cast the vertex coordinates to. - Defaults to the plain Python `float` used by the hardcoded - vertex coordinates. + :arg dtype: optional working dtype for tabulation. """ if spatial_dim == 0: cell = Point() @@ -1712,20 +1709,22 @@ def ufc_simplex(spatial_dim, dtype=None): else: raise RuntimeError(f"Can't create UFC simplex of dimension {spatial_dim}.") if dtype is not None: - cell.vertices = cast_vertices(cell.vertices, dtype) + cell.target_dtype = numpy.dtype(dtype) return cell def symmetric_simplex(spatial_dim, dtype=None): A = numpy.array([[2, 1, 1], [0, numpy.sqrt(3), numpy.sqrt(3)/3], - [0, 0, numpy.sqrt(6)*(2/3)]], dtype=dtype) + [0, 0, numpy.sqrt(6)*(2/3)]], dtype=float) A = A[:spatial_dim, :][:, :spatial_dim] b = A.sum(axis=1) * (-1 / (1 + spatial_dim)) Ref1 = ufc_simplex(spatial_dim, dtype=dtype) v = numpy.dot(Ref1.get_vertices(), A.T) + b[None, :] vertices = tuple(map(tuple, v)) - return SymmetricSimplex(Ref1.get_shape(), vertices, Ref1.get_topology()) + cell = SymmetricSimplex(Ref1.get_shape(), vertices, Ref1.get_topology()) + cell.target_dtype = numpy.dtype(float if dtype is None else dtype) + return cell def ufc_cell(cell, dtype=None): @@ -1739,7 +1738,9 @@ def ufc_cell(cell, dtype=None): if " * " in celltype: # Tensor product cell - return TensorProductCell(*(ufc_cell(c, dtype=dtype) for c in celltype.split(" * "))) + ref_el = TensorProductCell(*(ufc_cell(c, dtype=dtype) for c in celltype.split(" * "))) + ref_el.target_dtype = numpy.dtype(float if dtype is None else dtype) + return ref_el elif celltype == "quadrilateral": return ufc_hypercube(2, dtype=dtype) elif celltype == "hexahedron": diff --git a/test/FIAT/unit/test_hct.py b/test/FIAT/unit/test_hct.py index 73109f91..787c5849 100644 --- a/test/FIAT/unit/test_hct.py +++ b/test/FIAT/unit/test_hct.py @@ -10,8 +10,7 @@ @pytest.fixture(params=(numpy.float64, numpy.float32), ids=("float64", "float32")) def cell(request): K = ufc_simplex(2, dtype=request.param) - K.vertices = tuple(map(tuple, numpy.asarray( - ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84)), dtype=request.param))) + K.vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84)) return K diff --git a/test/FIAT/unit/test_powell_sabin.py b/test/FIAT/unit/test_powell_sabin.py index 31f19ed7..900adc1b 100644 --- a/test/FIAT/unit/test_powell_sabin.py +++ b/test/FIAT/unit/test_powell_sabin.py @@ -16,8 +16,8 @@ def cell(request): def test_powell_sabin_constant(cell, el): # Test that bfs associated with point evaluation sum up to 1 fe = el(cell) - assert (numpy.asarray(fe.get_reference_complex().vertices).dtype - == numpy.asarray(cell.vertices).dtype) + assert numpy.asarray(fe.get_reference_complex().vertices).dtype == numpy.float64 + assert fe.get_reference_complex().target_dtype == cell.target_dtype assert fe.get_coeffs().dtype == numpy.float64 pts = make_lattice(cell.get_vertices(), 3) diff --git a/test/finat/test_create_finat_element.py b/test/finat/test_create_finat_element.py index 66624713..4687a9b8 100644 --- a/test/finat/test_create_finat_element.py +++ b/test/finat/test_create_finat_element.py @@ -168,15 +168,17 @@ def test_cache_hit_vector(ufl_vector_element): def test_dtype_reaches_reference_cell(ufl_element): - """dtype passed to create_element must reach the constructed element's - reference cell, not silently fall back to FIAT's float64 default.""" + """Check that target dtype is independent of construction dtype.""" default = create_element(ufl_element) single = create_element(ufl_element, dtype=numpy.float32) double = create_element(ufl_element, dtype=numpy.float64) assert numpy.array(default.cell.vertices).dtype == numpy.float64 - assert numpy.array(single.cell.vertices).dtype == numpy.float32 + assert numpy.array(single.cell.vertices).dtype == numpy.float64 assert numpy.array(double.cell.vertices).dtype == numpy.float64 + assert default.cell.target_dtype == numpy.float64 + assert single.cell.target_dtype == numpy.float32 + assert double.cell.target_dtype == numpy.float64 def test_dtype_cache_distinguishes(ufl_element): From a9182aa0eec50fc8bf72117b93c16e712186979d Mon Sep 17 00:00:00 2001 From: Hardik Kothari Date: Fri, 4 Sep 2026 12:26:01 +0200 Subject: [PATCH 6/6] Update precision test for double construction --- test/FIAT/unit/test_precision.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/FIAT/unit/test_precision.py b/test/FIAT/unit/test_precision.py index 30b82da6..f797ecdd 100644 --- a/test/FIAT/unit/test_precision.py +++ b/test/FIAT/unit/test_precision.py @@ -46,7 +46,9 @@ def test_dtype_propagates_into_symbolic_tabulation(): fe64 = macro_element(dtype=numpy.float64) fe32 = macro_element(dtype=numpy.float32) assert numpy.array(fe64.ref_complex.vertices).dtype == numpy.float64 - assert numpy.array(fe32.ref_complex.vertices).dtype == numpy.float32 + assert numpy.array(fe32.ref_complex.vertices).dtype == numpy.float64 + assert fe64.ref_complex.target_dtype == numpy.float64 + assert fe32.ref_complex.target_dtype == numpy.float32 x0 = gem.Variable("x0", ()) tab64 = fe64.tabulate(0, (x0,))[(0,)][0]