From 998753ca88b849bd92b2b2b261b9cbdf2278e544 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 28 May 2024 18:44:52 +0100 Subject: [PATCH 1/2] Get interior nodes directly from recursivenodes --- FIAT/reference_element.py | 26 ++++++-------------------- 1 file changed, 6 insertions(+), 20 deletions(-) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 8ad5bd1a5..f7b22fa43 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -26,7 +26,7 @@ from math import factorial import numpy -from recursivenodes.nodes import _decode_family, _recursive +from recursivenodes import recursive_nodes from FIAT.orientation_utils import ( make_cell_orientation_reflection_map_simplex, @@ -41,20 +41,6 @@ TENSORPRODUCT = 99 -def multiindex_equal(d, isum, imin=0): - """A generator for d-tuple multi-indices whose sum is isum and minimum is imin. - """ - if d <= 0: - return - imax = isum - (d - 1) * imin - if imax < imin: - return - for i in range(imin, imax): - for a in multiindex_equal(d - 1, isum - i, imin=imin): - yield a + (i,) - yield (imin,) * (d - 1) + (imax,) - - def lattice_iter(start, finish, depth): """Generator iterating over the depth-dimensional lattice of integers between start and (finish-1). This works on simplices in @@ -85,11 +71,11 @@ def make_lattice(verts, n, interior=0, variant=None): "equispaced_interior": "equi_interior", "gll": "lgl"} family = recursivenodes_families.get(variant, variant) - family = _decode_family(family) - D = len(verts) - X = numpy.array(verts) - get_point = lambda alpha: tuple(numpy.dot(_recursive(D - 1, n, alpha, family), X)) - return list(map(get_point, multiindex_equal(D, n, interior))) + D = len(verts) - 1 + X = numpy.asarray(verts[::-1]) + bary = recursive_nodes(D, n, family=family, interior=interior) + pts = numpy.dot(bary, X) + return list(map(tuple, pts)) def linalg_subspace_intersection(A, B): From 15916e8ef381892a7e6d6c1d4d8f7d1b397ca8fb Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 28 May 2024 18:45:53 +0100 Subject: [PATCH 2/2] Fix some warnings --- FIAT/discontinuous_lagrange.py | 3 ++- FIAT/dual_set.py | 2 +- FIAT/hct.py | 2 +- FIAT/morley.py | 2 +- FIAT/orientation_utils.py | 3 ++- FIAT/reference_element.py | 2 +- test/unit/test_kong_mulder_veldhuizen.py | 3 +-- 7 files changed, 9 insertions(+), 8 deletions(-) diff --git a/FIAT/discontinuous_lagrange.py b/FIAT/discontinuous_lagrange.py index 85907b56b..cb2e932a0 100644 --- a/FIAT/discontinuous_lagrange.py +++ b/FIAT/discontinuous_lagrange.py @@ -7,6 +7,7 @@ import itertools import numpy as np +from math import factorial from FIAT import finite_element, polynomial_set, dual_set, functional, P0 from FIAT.reference_element import LINE, make_lattice from FIAT.orientation_utils import make_entity_permutations_simplex @@ -17,7 +18,7 @@ def make_entity_permutations(dim, npoints): if npoints <= 0: - return {o: [] for o in range(np.math.factorial(dim + 1))} + return {o: [] for o in range(factorial(dim + 1))} # DG nodes are numbered, in order of significance, # - by g0: entity dim (vertices first, then edges, then ...) # - by g1: entity ids (DoFs on entities of smaller ids first) diff --git a/FIAT/dual_set.py b/FIAT/dual_set.py index d7015996a..7360491c5 100644 --- a/FIAT/dual_set.py +++ b/FIAT/dual_set.py @@ -177,7 +177,7 @@ def to_riesz(self, poly_set): for pt, wac_list in ell.deriv_dict.items(): j = dpts.index(pt) for (w, alpha, c) in wac_list: - dwts[alpha][i][c][j] = w + dwts[alpha][(i, *c, j)] = w for alpha in dwts: mat[ells] += numpy.dot(dwts[alpha], dexpansion_values[alpha].T) return mat diff --git a/FIAT/hct.py b/FIAT/hct.py index 00668653d..d1b18952f 100644 --- a/FIAT/hct.py +++ b/FIAT/hct.py @@ -31,7 +31,7 @@ def __init__(self, ref_el, degree, reduced=False): k = 2 if reduced else 0 Q = create_quadrature(rline, degree-1+k) qpts = Q.get_points() - f_at_qpts = eval_jacobi(0, 0, k, 2.0*qpts - 1) + f_at_qpts = eval_jacobi(0, 0, k, 2.0*qpts[:, 0] - 1) for e in sorted(top[1]): cur = len(nodes) nodes.append(IntegralMomentOfNormalDerivative(ref_el, e, Q, f_at_qpts)) diff --git a/FIAT/morley.py b/FIAT/morley.py index 67896426d..bcaf142b1 100644 --- a/FIAT/morley.py +++ b/FIAT/morley.py @@ -41,7 +41,7 @@ def __init__(self, ref_el): degree = 2 Q = create_quadrature(rline, degree-1) qpts = Q.get_points() - scale = numpy.ones(qpts.shape) + scale = numpy.ones((len(qpts),)) entity_ids[1] = {} for e in sorted(top[1]): diff --git a/FIAT/orientation_utils.py b/FIAT/orientation_utils.py index 35ddc5293..5354df8c1 100644 --- a/FIAT/orientation_utils.py +++ b/FIAT/orientation_utils.py @@ -1,6 +1,7 @@ import itertools from collections.abc import Sequence import numpy as np +from math import factorial def make_entity_permutations_simplex(dim, npoints): @@ -53,7 +54,7 @@ def make_entity_permutations_simplex(dim, npoints): from FIAT.polynomial_set import mis if npoints <= 0: - return {o: [] for o in range(np.math.factorial(dim + 1))} + return {o: [] for o in range(factorial(dim + 1))} a = np.array(sorted(mis(dim + 1, npoints - 1)), dtype=int)[:, ::-1] index_perms = sorted(itertools.permutations(range(dim + 1))) perms = {} diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index f7b22fa43..f941a1c79 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -546,7 +546,7 @@ def is_simplex(self): return True def symmetry_group_size(self, dim): - return numpy.math.factorial(dim + 1) + return factorial(dim + 1) def cell_orientation_reflection_map(self): """Return the map indicating whether each possible cell orientation causes reflection (``1``) or not (``0``).""" diff --git a/test/unit/test_kong_mulder_veldhuizen.py b/test/unit/test_kong_mulder_veldhuizen.py index 5fce238ec..40d09e7bf 100644 --- a/test/unit/test_kong_mulder_veldhuizen.py +++ b/test/unit/test_kong_mulder_veldhuizen.py @@ -1,5 +1,6 @@ import numpy as np import pytest +from math import factorial as fct from FIAT.reference_element import UFCInterval, UFCTriangle, UFCTetrahedron from FIAT import create_quadrature, make_quadrature, polynomial_set @@ -12,7 +13,6 @@ @pytest.mark.parametrize("p_d", [(1, 1), (2, 3), (3, 4)]) def test_kmv_quad_tet_schemes(p_d): # noqa: W503 - fct = np.math.factorial p, d = p_d q = create_quadrature(Te, p, "KMV") for i in range(d + 1): @@ -30,7 +30,6 @@ def test_kmv_quad_tet_schemes(p_d): # noqa: W503 @pytest.mark.parametrize("p_d", [(1, 1), (2, 3), (3, 5), (4, 7), (5, 9)]) def test_kmv_quad_tri_schemes(p_d): - fct = np.math.factorial p, d = p_d q = create_quadrature(T, p, "KMV") for i in range(d + 1):