From 81e8177e13c0a0cf609a84af5459814f5fcceb73 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 10 Jul 2024 17:18:44 +0100 Subject: [PATCH 1/3] GLL: ensure natural ordering of the DOFs --- FIAT/gauss_lobatto_legendre.py | 2 +- FIAT/lagrange.py | 55 ++++++++++++++++-------- test/unit/test_fiat.py | 4 +- test/unit/test_gauss_lobatto_legendre.py | 3 +- 4 files changed, 41 insertions(+), 23 deletions(-) diff --git a/FIAT/gauss_lobatto_legendre.py b/FIAT/gauss_lobatto_legendre.py index d4589da0e..64aee5478 100644 --- a/FIAT/gauss_lobatto_legendre.py +++ b/FIAT/gauss_lobatto_legendre.py @@ -14,4 +14,4 @@ class GaussLobattoLegendre(lagrange.Lagrange): """Simplicial continuous element with nodes at the (recursive) Gauss-Lobatto-Legendre points.""" def __init__(self, ref_el, degree): - super(GaussLobattoLegendre, self).__init__(ref_el, degree, variant="gll") + super(GaussLobattoLegendre, self).__init__(ref_el, degree, variant="gll", lexsort=True) diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index 1201369c9..e04055bac 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -8,8 +8,9 @@ from FIAT import finite_element, polynomial_set, dual_set, functional from FIAT.orientation_utils import make_entity_permutations_simplex from FIAT.barycentric_interpolation import LagrangePolynomialSet, get_lagrange_points -from FIAT.reference_element import LINE +from FIAT.reference_element import multiindex_equal, make_lattice, LINE from FIAT.check_format_variant import parse_lagrange_variant +import numpy class LagrangeDualSet(dual_set.DualSet): @@ -17,28 +18,45 @@ class LagrangeDualSet(dual_set.DualSet): simplices of any dimension. Nodes are point evaluation at recursively-defined points. """ - def __init__(self, ref_el, degree, point_variant="equispaced"): + def __init__(self, ref_el, degree, point_variant="equispaced", lexsort=False): entity_ids = {} - nodes = [] - entity_permutations = {} - - # make nodes by getting points - # need to do this dimension-by-dimension, facet-by-facet top = ref_el.get_topology() + if lexsort: + # make nodes by getting points on a lattice + points = make_lattice(ref_el.get_vertices(), degree, variant=point_variant) + nodes = [functional.PointEvaluation(ref_el, x) for x in points] + + # map lattice coordinates to node indices + sd = ref_el.get_spatial_dimension() + alphas = multiindex_equal(sd+1, degree) + ids = dict(zip(alphas, range(len(points)))) - cur = 0 + # associate entities to node indices + ref_vertices = [tuple(int(i == j) for j in range(sd+1)) for i in range(sd+1)] + for dim in sorted(top): + entity_ids[dim] = {} + base_alphas = list(multiindex_equal(dim+1, degree, imin=1)) + for entity in sorted(top[dim]): + verts = [ref_vertices[i] for i in top[dim][entity]] + alphas = [] if len(base_alphas) == 0 else numpy.dot(base_alphas, verts) + entity_ids[dim][entity] = [ids[tuple(alpha)] for alpha in alphas] + else: + # make nodes by getting points + # need to do this dimension-by-dimension, facet-by-facet + nodes = [] + for dim in sorted(top): + entity_ids[dim] = {} + for entity in sorted(top[dim]): + cur = len(nodes) + pts_cur = ref_el.make_points(dim, entity, degree, variant=point_variant) + nodes.extend(functional.PointEvaluation(ref_el, x) for x in pts_cur) + entity_ids[dim][entity] = list(range(cur, len(nodes))) + + entity_permutations = {} for dim in sorted(top): - entity_ids[dim] = {} entity_permutations[dim] = {} perms = {0: [0]} if dim == 0 else make_entity_permutations_simplex(dim, degree - dim) for entity in sorted(top[dim]): - pts_cur = ref_el.make_points(dim, entity, degree, variant=point_variant) - nodes_cur = [functional.PointEvaluation(ref_el, x) - for x in pts_cur] - nnodes_cur = len(nodes_cur) - nodes += nodes_cur - entity_ids[dim][entity] = list(range(cur, cur + nnodes_cur)) - cur += nnodes_cur entity_permutations[dim][entity] = perms super(LagrangeDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) @@ -58,12 +76,13 @@ class Lagrange(finite_element.CiarletElement): variant='equispaced,Iso(2)' with degree=1 gives the P2:P1 iso element. variant='Alfeld' can be used to obtain a barycentrically refined macroelement for Scott-Vogelius. + :arg lexsort: Are we sorting the nodes lexicographically? """ - def __init__(self, ref_el, degree, variant="equispaced"): + def __init__(self, ref_el, degree, variant="equispaced", lexsort=False): splitting, point_variant = parse_lagrange_variant(variant) if splitting is not None: ref_el = splitting(ref_el) - dual = LagrangeDualSet(ref_el, degree, point_variant=point_variant) + dual = LagrangeDualSet(ref_el, degree, point_variant=point_variant, lexsort=lexsort) if ref_el.shape == LINE: # In 1D we can use the primal basis as the expansion set, # avoiding any round-off coming from a basis transformation diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index f81b8aa3f..65abfaef9 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -393,8 +393,8 @@ def test_empty_bubble(): @pytest.mark.parametrize('elements', [ (Lagrange(I, 2), Lagrange(I, 1), Bubble(I, 2)), - (GaussLobattoLegendre(I, 3), Lagrange(I, 1), - RestrictedElement(GaussLobattoLegendre(I, 3), restriction_domain="interior")), + (Lagrange(I, 3, variant="gll"), Lagrange(I, 1), + RestrictedElement(Lagrange(I, 3, variant="gll"), restriction_domain="interior")), (RaviartThomas(T, 2), RestrictedElement(RaviartThomas(T, 2), restriction_domain='facet'), RestrictedElement(RaviartThomas(T, 2), restriction_domain='interior')), diff --git a/test/unit/test_gauss_lobatto_legendre.py b/test/unit/test_gauss_lobatto_legendre.py index 37abe5a13..4c288d654 100644 --- a/test/unit/test_gauss_lobatto_legendre.py +++ b/test/unit/test_gauss_lobatto_legendre.py @@ -59,8 +59,7 @@ def test_edge_dofs(dim, degree): # Test that edge DOFs are located at the GLL quadrature points line = s if dim == 1 else s.construct_subelement(1) lr = quadrature.GaussLobattoLegendreQuadratureLineRule(line, degree + 1) - # Edge DOFs are ordered with the two vertex DOFs followed by the interior DOFs - quadrature_points = lr.pts[::degree] + lr.pts[1:-1] + quadrature_points = lr.get_points() entity_dofs = fe.entity_closure_dofs() edge_dofs = entity_dofs[1] From a9e96d951889d2e86150cc807dcfeffccefec53c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 17 Jul 2024 12:19:11 +0100 Subject: [PATCH 2/3] Lagrange: Always sort DOFs by entity. Entities may be sorted by dimension or by their support vertex ids. --- FIAT/gauss_lobatto_legendre.py | 2 +- FIAT/lagrange.py | 72 ++++++++++++++++------------------ 2 files changed, 34 insertions(+), 40 deletions(-) diff --git a/FIAT/gauss_lobatto_legendre.py b/FIAT/gauss_lobatto_legendre.py index 64aee5478..1627d3f09 100644 --- a/FIAT/gauss_lobatto_legendre.py +++ b/FIAT/gauss_lobatto_legendre.py @@ -14,4 +14,4 @@ class GaussLobattoLegendre(lagrange.Lagrange): """Simplicial continuous element with nodes at the (recursive) Gauss-Lobatto-Legendre points.""" def __init__(self, ref_el, degree): - super(GaussLobattoLegendre, self).__init__(ref_el, degree, variant="gll", lexsort=True) + super(GaussLobattoLegendre, self).__init__(ref_el, degree, variant="gll", lexsort_entities=True) diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index e04055bac..6ec47ccb3 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -8,57 +8,48 @@ from FIAT import finite_element, polynomial_set, dual_set, functional from FIAT.orientation_utils import make_entity_permutations_simplex from FIAT.barycentric_interpolation import LagrangePolynomialSet, get_lagrange_points -from FIAT.reference_element import multiindex_equal, make_lattice, LINE +from FIAT.reference_element import LINE from FIAT.check_format_variant import parse_lagrange_variant -import numpy class LagrangeDualSet(dual_set.DualSet): """The dual basis for Lagrange elements. This class works for - simplices of any dimension. Nodes are point evaluation at + simplicial complexes of any dimension. Nodes are point evaluation at recursively-defined points. + + :arg ref_el: The simplicial complex. + :arg degree: The polynomial degree. + :point_variant: The point distribution variant passed on to recursivenodes. + :lexsort_entities: A flag to sort entities by lexicographically by coordinate, + if false then entities are sorted first by dimension and then by + entity id. The DOFs are always sorted by the entity ordering + and then lexicographically by coordinate. """ - def __init__(self, ref_el, degree, point_variant="equispaced", lexsort=False): + def __init__(self, ref_el, degree, point_variant="equispaced", lexsort_entities=False): + nodes = [] entity_ids = {} - top = ref_el.get_topology() - if lexsort: - # make nodes by getting points on a lattice - points = make_lattice(ref_el.get_vertices(), degree, variant=point_variant) - nodes = [functional.PointEvaluation(ref_el, x) for x in points] - - # map lattice coordinates to node indices - sd = ref_el.get_spatial_dimension() - alphas = multiindex_equal(sd+1, degree) - ids = dict(zip(alphas, range(len(points)))) - - # associate entities to node indices - ref_vertices = [tuple(int(i == j) for j in range(sd+1)) for i in range(sd+1)] - for dim in sorted(top): - entity_ids[dim] = {} - base_alphas = list(multiindex_equal(dim+1, degree, imin=1)) - for entity in sorted(top[dim]): - verts = [ref_vertices[i] for i in top[dim][entity]] - alphas = [] if len(base_alphas) == 0 else numpy.dot(base_alphas, verts) - entity_ids[dim][entity] = [ids[tuple(alpha)] for alpha in alphas] - else: - # make nodes by getting points - # need to do this dimension-by-dimension, facet-by-facet - nodes = [] - for dim in sorted(top): - entity_ids[dim] = {} - for entity in sorted(top[dim]): - cur = len(nodes) - pts_cur = ref_el.make_points(dim, entity, degree, variant=point_variant) - nodes.extend(functional.PointEvaluation(ref_el, x) for x in pts_cur) - entity_ids[dim][entity] = list(range(cur, len(nodes))) - entity_permutations = {} + top = ref_el.get_topology() for dim in sorted(top): + entity_ids[dim] = {} entity_permutations[dim] = {} perms = {0: [0]} if dim == 0 else make_entity_permutations_simplex(dim, degree - dim) for entity in sorted(top[dim]): entity_permutations[dim][entity] = perms + entities = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim])] + if lexsort_entities: + # sort the entities by their support vertex ids + support = [top[dim][entity] for dim, entity in entities] + entities = [entity for verts, entity in sorted(zip(support, entities))] + + # make nodes by getting points + # need to do this entity-by-entity + for dim, entity in entities: + cur = len(nodes) + pts_cur = ref_el.make_points(dim, entity, degree, variant=point_variant) + nodes.extend(functional.PointEvaluation(ref_el, x) for x in pts_cur) + entity_ids[dim][entity] = list(range(cur, len(nodes))) super(LagrangeDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations) @@ -76,13 +67,16 @@ class Lagrange(finite_element.CiarletElement): variant='equispaced,Iso(2)' with degree=1 gives the P2:P1 iso element. variant='Alfeld' can be used to obtain a barycentrically refined macroelement for Scott-Vogelius. - :arg lexsort: Are we sorting the nodes lexicographically? + :arg lexsort_entities: A flag to sort entities lexicographically by coordinate, + if false then entities are sorted first by dimension and then by + entity id. The DOFs are always sorted by the entity ordering + and then lexicographically by coordinate. """ - def __init__(self, ref_el, degree, variant="equispaced", lexsort=False): + def __init__(self, ref_el, degree, variant="equispaced", lexsort_entities=False): splitting, point_variant = parse_lagrange_variant(variant) if splitting is not None: ref_el = splitting(ref_el) - dual = LagrangeDualSet(ref_el, degree, point_variant=point_variant, lexsort=lexsort) + dual = LagrangeDualSet(ref_el, degree, point_variant=point_variant, lexsort_entities=lexsort_entities) if ref_el.shape == LINE: # In 1D we can use the primal basis as the expansion set, # avoiding any round-off coming from a basis transformation From 3d6db424dcdf4ab6c7affccb1c54fd401e5ba760 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 17 Jul 2024 15:42:32 +0100 Subject: [PATCH 3/3] Fix docs --- FIAT/gauss_lobatto_legendre.py | 2 +- FIAT/lagrange.py | 28 ++++++++++++++-------------- 2 files changed, 15 insertions(+), 15 deletions(-) diff --git a/FIAT/gauss_lobatto_legendre.py b/FIAT/gauss_lobatto_legendre.py index 1627d3f09..e97c7d886 100644 --- a/FIAT/gauss_lobatto_legendre.py +++ b/FIAT/gauss_lobatto_legendre.py @@ -14,4 +14,4 @@ class GaussLobattoLegendre(lagrange.Lagrange): """Simplicial continuous element with nodes at the (recursive) Gauss-Lobatto-Legendre points.""" def __init__(self, ref_el, degree): - super(GaussLobattoLegendre, self).__init__(ref_el, degree, variant="gll", lexsort_entities=True) + super(GaussLobattoLegendre, self).__init__(ref_el, degree, variant="gll", sort_entities=True) diff --git a/FIAT/lagrange.py b/FIAT/lagrange.py index 6ec47ccb3..51a996154 100644 --- a/FIAT/lagrange.py +++ b/FIAT/lagrange.py @@ -19,13 +19,13 @@ class LagrangeDualSet(dual_set.DualSet): :arg ref_el: The simplicial complex. :arg degree: The polynomial degree. - :point_variant: The point distribution variant passed on to recursivenodes. - :lexsort_entities: A flag to sort entities by lexicographically by coordinate, - if false then entities are sorted first by dimension and then by - entity id. The DOFs are always sorted by the entity ordering - and then lexicographically by coordinate. + :arg point_variant: The point distribution variant passed on to recursivenodes. + :arg sort_entities: A flag to sort entities by support vertex ids. + If false then entities are sorted first by dimension and then by + entity id. The DOFs are always sorted by the entity ordering + and then lexicographically by lattice multiindex. """ - def __init__(self, ref_el, degree, point_variant="equispaced", lexsort_entities=False): + def __init__(self, ref_el, degree, point_variant="equispaced", sort_entities=False): nodes = [] entity_ids = {} entity_permutations = {} @@ -38,8 +38,8 @@ def __init__(self, ref_el, degree, point_variant="equispaced", lexsort_entities= entity_permutations[dim][entity] = perms entities = [(dim, entity) for dim in sorted(top) for entity in sorted(top[dim])] - if lexsort_entities: - # sort the entities by their support vertex ids + if sort_entities: + # sort the entities by support vertex ids support = [top[dim][entity] for dim, entity in entities] entities = [entity for verts, entity in sorted(zip(support, entities))] @@ -67,16 +67,16 @@ class Lagrange(finite_element.CiarletElement): variant='equispaced,Iso(2)' with degree=1 gives the P2:P1 iso element. variant='Alfeld' can be used to obtain a barycentrically refined macroelement for Scott-Vogelius. - :arg lexsort_entities: A flag to sort entities lexicographically by coordinate, - if false then entities are sorted first by dimension and then by - entity id. The DOFs are always sorted by the entity ordering - and then lexicographically by coordinate. + :arg sort_entities: A flag to sort entities by support vertex ids. + If false then entities are sorted first by dimension and then by + entity id. The DOFs are always sorted by the entity ordering + and then lexicographically by lattice multiindex. """ - def __init__(self, ref_el, degree, variant="equispaced", lexsort_entities=False): + def __init__(self, ref_el, degree, variant="equispaced", sort_entities=False): splitting, point_variant = parse_lagrange_variant(variant) if splitting is not None: ref_el = splitting(ref_el) - dual = LagrangeDualSet(ref_el, degree, point_variant=point_variant, lexsort_entities=lexsort_entities) + dual = LagrangeDualSet(ref_el, degree, point_variant=point_variant, sort_entities=sort_entities) if ref_el.shape == LINE: # In 1D we can use the primal basis as the expansion set, # avoiding any round-off coming from a basis transformation