From 69073c8f6bd4217c3a8800fa01489dcd9d7d39f8 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 26 Mar 2026 12:08:12 +0000 Subject: [PATCH 1/3] Hermite on 1D manifolds --- FIAT/expansions.py | 2 +- FIAT/hermite.py | 57 +++++++++++---------------------- FIAT/reference_element.py | 5 +++ finat/hermite.py | 24 ++++++++------ test/finat/conftest.py | 5 +-- test/finat/test_zany_mapping.py | 8 +++++ 6 files changed, 51 insertions(+), 50 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 0ddea371a..1473d9d6a 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -580,7 +580,7 @@ def __init__(self, ref_el, **kwargs): def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): """Returns a dict of tabulations such that tabulations[alpha][i, j] = D^alpha phi_i(pts[j]).""" - if self.variant is not None: + if self.variant is not None or self.ref_el.get_spatial_dimension() != 1: return super()._tabulate_on_cell(n, pts, order=order, cell=cell, direction=direction) A, b = self.affine_mappings[cell] diff --git a/FIAT/hermite.py b/FIAT/hermite.py index 79c83ff97..d1342d3e7 100644 --- a/FIAT/hermite.py +++ b/FIAT/hermite.py @@ -9,60 +9,41 @@ class CubicHermiteDualSet(dual_set.DualSet): - """The dual basis for Lagrange elements. This class works for - simplices of any dimension. Nodes are point evaluation at - equispaced points.""" + """The dual basis for Hermite elements. This class works for + simplices of any dimension. Nodes are first order jet at + vertices and point evaluation at barycenters of 2D entities.""" def __init__(self, ref_el): - entity_ids = {} - nodes = [] - cur = 0 - # make nodes by getting points # need to do this dimension-by-dimension, facet-by-facet top = ref_el.get_topology() verts = ref_el.get_vertices() - sd = ref_el.get_spatial_dimension() - - # get jet at each vertex + sd = ref_el.get_topological_dimension() + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] - entity_ids[0] = {} + # get first order jet at each vertex for v in sorted(top[0]): + cur = len(nodes) nodes.append(functional.PointEvaluation(ref_el, verts[v])) - pd = functional.PointDerivative - for i in range(sd): - alpha = [0] * sd - alpha[i] = 1 - - nodes.append(pd(ref_el, verts[v], alpha)) + if sd == 1: + # in 1D use normal derivative to support manifolds + nodes.append(functional.PointNormalDerivative(ref_el, v, verts[v])) + else: + nodes.extend(functional.PointDerivative(ref_el, verts[v], alpha) + for alpha in polynomial_set.mis(sd, 1)) + entity_ids[0][v].extend(range(cur, len(nodes))) - entity_ids[0][v] = list(range(cur, cur + 1 + sd)) - cur += sd + 1 - - # now only have dofs at the barycenter, which is the - # maximal dimension # no edge dof - - entity_ids[1] = {} - for i in top[1]: - entity_ids - entity_ids[1][i] = [] - + # now only add dofs at the barycenter of the 2D entities if sd > 1: # face dof # point evaluation at barycenter - entity_ids[2] = {} for f in sorted(top[2]): + cur = len(nodes) pt = ref_el.make_points(2, f, 3)[0] - n = functional.PointEvaluation(ref_el, pt) - nodes.append(n) - entity_ids[2][f] = list(range(cur, cur + 1)) - cur += 1 - - for dim in range(3, sd + 1): - entity_ids[dim] = {} - for facet in top[dim]: - entity_ids[dim][facet] = [] + nodes.append(functional.PointEvaluation(ref_el, pt)) + entity_ids[2][f].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index b896059e1..becf0d7c2 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -1001,6 +1001,11 @@ def __init__(self): 1: edges} super().__init__(LINE, verts, topology) + def compute_normal(self, i): + "UFC consistent normal" + n = self.compute_tangents(1, 0)[0] + return n / numpy.linalg.norm(n) + class DefaultTriangle(DefaultSimplex): """This is the reference triangle with vertices (-1.0,-1.0), diff --git a/finat/hermite.py b/finat/hermite.py index 5b3a28115..cfce7f917 100644 --- a/finat/hermite.py +++ b/finat/hermite.py @@ -1,5 +1,6 @@ import FIAT -from gem import ListTensor +import numpy +from gem import ListTensor, partial_indexed from finat.citations import cite from finat.fiat_elements import ScalarFiatElement @@ -15,18 +16,23 @@ def basis_transformation(self, coordinate_mapping): Js = [coordinate_mapping.jacobian_at(vertex) for vertex in self.cell.get_vertices()] + pns = coordinate_mapping.physical_normals() h = coordinate_mapping.cell_size() - d = self.cell.get_dimension() M = identity(self.space_dimension()) - cur = 0 - for i in range(d+1): - cur += 1 # skip the vertex + entity_ids = self.entity_dofs() + for i in entity_ids[0]: + # skip the PointEvaluation DOF + vids = entity_ids[0][i][1:] J = Js[i] - for j in range(d): - for k in range(d): - M[cur+j, cur+k] = J[j, k] / h[i] - cur += d + + gdim, tdim = J.shape + if gdim != tdim: + assert tdim == 1 + J = partial_indexed(pns, (i,)) @ J + + Jnp = numpy.reshape([J[i] for i in numpy.ndindex(J.shape)], J.shape) + M[numpy.ix_(vids, vids)] = Jnp * (1 / h[i]) return ListTensor(M) diff --git a/test/finat/conftest.py b/test/finat/conftest.py index 1fab8dd3d..a4d875253 100644 --- a/test/finat/conftest.py +++ b/test/finat/conftest.py @@ -107,18 +107,19 @@ def scaled_simplex(dim, scale): @pytest.fixture def ref_el(): - K = {dim: FIAT.ufc_simplex(dim) for dim in (2, 3)} + K = {dim: FIAT.ufc_simplex(dim) for dim in (1, 2, 3)} return K @pytest.fixture def phys_el(): - K = {dim: FIAT.ufc_simplex(int(dim)) for dim in (2, 2.5, 3)} + K = {dim: FIAT.ufc_simplex(int(dim)) for dim in (1.5, 2, 2.5, 3)} K[2].vertices = ((0.0, 0.1), (1.17, -0.09), (0.15, 1.84)) K[3].vertices = ((0, 0, 0), (1., 0.1, -0.37), (0.01, 0.987, -.23), (-0.1, -0.2, 1.38)) + K[1.5].vertices = K[2].vertices[:-1] K[2.5].vertices = K[3].vertices[:-1] return K diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index fe74d2ad4..92e20717d 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -119,6 +119,13 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), pp.pformat((np.round(error, 8).tolist(), *inds)) +@pytest.mark.parametrize("element", [ + finat.Hermite, + ]) +def test_C1_interval(ref_to_phys, element): + check_zany_mapping(element, ref_to_phys[1.5]) + + @pytest.mark.parametrize("element", [ finat.Morley, finat.Hermite, @@ -132,6 +139,7 @@ def test_C1_triangle(ref_to_phys, element): @pytest.mark.parametrize("element", [ finat.Morley, + finat.Hermite, finat.Walkington, ]) def test_C1_tetrahedron(ref_to_phys, element): From ad1f782e2520f990950d33d9ea4ffefd55bb3206 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 26 Mar 2026 14:48:04 +0000 Subject: [PATCH 2/3] Higher-order Hermite in 1D --- FIAT/hermite.py | 41 +++++++++++++++++++-------------- finat/hermite.py | 4 ++-- test/finat/test_zany_mapping.py | 10 ++++---- 3 files changed, 31 insertions(+), 24 deletions(-) diff --git a/FIAT/hermite.py b/FIAT/hermite.py index d1342d3e7..970887f3e 100644 --- a/FIAT/hermite.py +++ b/FIAT/hermite.py @@ -13,35 +13,41 @@ class CubicHermiteDualSet(dual_set.DualSet): simplices of any dimension. Nodes are first order jet at vertices and point evaluation at barycenters of 2D entities.""" - def __init__(self, ref_el): + def __init__(self, ref_el, degree, variant=None): # make nodes by getting points # need to do this dimension-by-dimension, facet-by-facet top = ref_el.get_topology() - verts = ref_el.get_vertices() sd = ref_el.get_topological_dimension() entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} nodes = [] # get first order jet at each vertex for v in sorted(top[0]): + pt, = ref_el.make_points(0, v, degree, variant=variant) cur = len(nodes) - nodes.append(functional.PointEvaluation(ref_el, verts[v])) + nodes.append(functional.PointEvaluation(ref_el, pt)) if sd == 1: - # in 1D use normal derivative to support manifolds - nodes.append(functional.PointNormalDerivative(ref_el, v, verts[v])) + # use normal derivative to support manifolds in 1D + nodes.append(functional.PointNormalDerivative(ref_el, v, pt)) else: - nodes.extend(functional.PointDerivative(ref_el, verts[v], alpha) + nodes.extend(functional.PointDerivative(ref_el, pt, alpha) for alpha in polynomial_set.mis(sd, 1)) entity_ids[0][v].extend(range(cur, len(nodes))) - # no edge dof - # now only add dofs at the barycenter of the 2D entities - if sd > 1: - # face dof - # point evaluation at barycenter + if sd == 1: + # edge dofs: point evaluations to support higher order in 1D + for e in sorted(top[1]): + cur = len(nodes) + pts = ref_el.make_points(1, e, degree-2, variant=variant) + nodes.extend(functional.PointEvaluation(ref_el, pt) for pt in pts) + entity_ids[1][e].extend(range(cur, len(nodes))) + else: + assert degree == 3 + # no edge dof + # face dof: point evaluation at barycenter for f in sorted(top[2]): cur = len(nodes) - pt = ref_el.make_points(2, f, 3)[0] + pt, = ref_el.make_points(2, f, degree, variant=variant) nodes.append(functional.PointEvaluation(ref_el, pt)) entity_ids[2][f].extend(range(cur, len(nodes))) @@ -51,9 +57,10 @@ def __init__(self, ref_el): class CubicHermite(finite_element.CiarletElement): """The cubic Hermite finite element. It is what it is.""" - def __init__(self, ref_el, deg=3): - assert deg == 3 - poly_set = polynomial_set.ONPolynomialSet(ref_el, 3) - dual = CubicHermiteDualSet(ref_el) + def __init__(self, ref_el, degree=3, variant=None): + if variant is None: + variant = "gll" + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) + dual = CubicHermiteDualSet(ref_el, degree) - super().__init__(poly_set, dual, 3) + super().__init__(poly_set, dual, degree) diff --git a/finat/hermite.py b/finat/hermite.py index cfce7f917..61262c0ab 100644 --- a/finat/hermite.py +++ b/finat/hermite.py @@ -8,9 +8,9 @@ class Hermite(PhysicallyMappedElement, ScalarFiatElement): - def __init__(self, cell, degree=3): + def __init__(self, cell, degree=3, variant=None): cite("Ciarlet1972") - super().__init__(FIAT.CubicHermite(cell)) + super().__init__(FIAT.CubicHermite(cell, degree=degree, variant=variant)) def basis_transformation(self, coordinate_mapping): Js = [coordinate_mapping.jacobian_at(vertex) diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 92e20717d..eb22eb82a 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -119,11 +119,11 @@ def check_zany_mapping(element, ref_to_phys, *args, **kwargs): assert np.allclose(ref_vals_zany, phys_vals[:num_dofs]), pp.pformat((np.round(error, 8).tolist(), *inds)) -@pytest.mark.parametrize("element", [ - finat.Hermite, - ]) -def test_C1_interval(ref_to_phys, element): - check_zany_mapping(element, ref_to_phys[1.5]) +@pytest.mark.parametrize("element, degree", [ + *((finat.Hermite, k) for k in range(3, 6)), +]) +def test_C1_interval(ref_to_phys, element, degree): + check_zany_mapping(element, ref_to_phys[1.5], degree) @pytest.mark.parametrize("element", [ From 05aabd6a19ca1743b93087fa6d665ac8095d3bbd Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 6 May 2026 15:17:42 +0100 Subject: [PATCH 3/3] pass variant --- FIAT/hermite.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FIAT/hermite.py b/FIAT/hermite.py index 970887f3e..15bb4eeee 100644 --- a/FIAT/hermite.py +++ b/FIAT/hermite.py @@ -61,6 +61,6 @@ def __init__(self, ref_el, degree=3, variant=None): if variant is None: variant = "gll" poly_set = polynomial_set.ONPolynomialSet(ref_el, degree) - dual = CubicHermiteDualSet(ref_el, degree) + dual = CubicHermiteDualSet(ref_el, degree, variant=variant) super().__init__(poly_set, dual, degree)