Skip to content

Commit

Permalink
revert some changes
Browse files Browse the repository at this point in the history
  • Loading branch information
@pbrubeck
pbrubeck committed Nov 6, 2023
1 parent dc98d15 commit 79a47d0
Show file tree
Hide file tree
Showing 6 changed files with 44 additions and 47 deletions.
4 changes: 1 addition & 3 deletions FIAT/expansions.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,9 +75,7 @@ def dubiner_recurrence(dim, n, ref_pts, phi, jacobian=None, dphi=None):
for codim in range(dim):
# Extend the basis from codim to codim + 1
fa, fb, fc, dfa, dfb, dfc = jacobi_factors(*X[codim:codim+3], *dX[codim:codim+3])
# Get indices of low-dimensional basis
alphas = [tuple()] if codim == 0 else reference_element.lattice_iter(0, n, codim)
for sub_index in alphas:
for sub_index in reference_element.lattice_iter(0, n, codim):
# handle i = 1
icur = idx(*sub_index, 0)
inext = idx(*sub_index, 1)
Expand Down
2 changes: 1 addition & 1 deletion FIAT/gauss_legendre.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -32,7 +32,7 @@ def __init__(self, ref_el, degree):
entity_permutations[dim][entity] = perms

# make nodes by getting points
pts = make_lattice(ref_el.get_vertices(), degree, family="gl")
pts = make_lattice(ref_el.get_vertices(), degree, variant="gl")
nodes = [functional.PointEvaluation(ref_el, x) for x in pts]
entity_ids[dim][0] = list(range(len(nodes)))
super(GaussLegendreDualSet, self).__init__(nodes, ref_el, entity_ids, entity_permutations)
Expand Down
2 changes: 1 addition & 1 deletion FIAT/gauss_lobatto_legendre.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@ class GaussLobattoLegendre(finite_element.CiarletElement):
def __init__(self, ref_el, degree):
if ref_el.shape not in {LINE, TRIANGLE, TETRAHEDRON}:
raise ValueError("Gauss-Lobatto-Legendre elements are only defined on simplices.")
dual = lagrange.LagrangeDualSet(ref_el, degree, family="lgl")
dual = lagrange.LagrangeDualSet(ref_el, degree, variant="gll")
if ref_el.shape == LINE:
points = []
for node in dual.nodes:
Expand Down
4 changes: 2 additions & 2 deletions FIAT/lagrange.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@ class LagrangeDualSet(dual_set.DualSet):
simplices of any dimension. Nodes are point evaluation at
equispaced points."""

def __init__(self, ref_el, degree, family=None):
def __init__(self, ref_el, degree, variant=None):
entity_ids = {}
nodes = []
entity_permutations = {}
Expand All @@ -29,7 +29,7 @@ def __init__(self, ref_el, degree, family=None):
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, family=family)
pts_cur = ref_el.make_points(dim, entity, degree, variant=variant)
nodes_cur = [functional.PointEvaluation(ref_el, x)
for x in pts_cur]
nnodes_cur = len(nodes_cur)
Expand Down
25 changes: 12 additions & 13 deletions FIAT/reference_element.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,29 +57,28 @@ def multiindex_equal(d, isum, imin=0):
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
1d, 2d, 3d, and beyond"""
0d, 1d, 2d, 3d, and beyond"""
if depth == 0:
return
elif depth == 1:
for ii in range(start, finish):
yield [ii]
yield tuple()
else:
for ii in range(start, finish):
for jj in lattice_iter(start, finish - ii, depth - 1):
yield jj + [ii]
yield jj + (ii,)


def make_lattice(verts, n, interior=0, family=None):
def make_lattice(verts, n, interior=0, variant=None):
"""Constructs a lattice of points on the simplex defined by verts.
For example, the 1:st order lattice will be just the vertices.
The optional argument interior specifies how many points from
the boundary to omit. For example, on a line with n = 2,
and interior = 0, this function will return the vertices and
midpoint, but with interior = 1, it will only return the
midpoint."""
if family is None or family == "equispaced":
family = "equi"
family = _decode_family(family)
if variant is None or variant == "equispaced":
variant = "equi"
elif variant == "gll":
variant = "lgl"
family = _decode_family(variant)
D = len(verts)
X = numpy.array(verts)
get_point = lambda alpha: tuple(numpy.dot(_recursive(D - 1, n, alpha, family), X))
Expand Down Expand Up @@ -406,7 +405,7 @@ def compute_face_edge_tangents(self, dim, entity_id):
edge_ts.append(vert_coords[dest] - vert_coords[source])
return edge_ts

def make_points(self, dim, entity_id, order, family=None):
def make_points(self, dim, entity_id, order, variant=None):
"""Constructs a lattice of points on the entity_id:th
facet of dimension dim. Order indicates how many points to
include in each direction."""
Expand All @@ -416,9 +415,9 @@ def make_points(self, dim, entity_id, order, family=None):
entity_verts = \
self.get_vertices_of_subcomplex(
self.get_topology()[dim][entity_id])
return make_lattice(entity_verts, order, 1, family=family)
return make_lattice(entity_verts, order, 1, variant=variant)
elif dim == self.get_spatial_dimension():
return make_lattice(self.get_vertices(), order, 1, family=family)
return make_lattice(self.get_vertices(), order, 1, variant=variant)
else:
raise ValueError("illegal dimension")

Expand Down
54 changes: 27 additions & 27 deletions test/unit/test_fiat.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -150,12 +150,12 @@ def __init__(self, a, b):
'RaviartThomas(S, 1, variant="integral(2)")',
'RaviartThomas(S, 2, variant="integral(3)")',
'RaviartThomas(S, 3, variant="integral(4)")',
'RaviartThomas(T, 1, variant="point")',
'RaviartThomas(T, 2, variant="point")',
'RaviartThomas(T, 3, variant="point")',
'RaviartThomas(S, 1, variant="point")',
'RaviartThomas(S, 2, variant="point")',
'RaviartThomas(S, 3, variant="point")',
'RaviartThomas(T, 1, variant="equispaced")',
'RaviartThomas(T, 2, variant="equispaced")',
'RaviartThomas(T, 3, variant="equispaced")',
'RaviartThomas(S, 1, variant="equispaced")',
'RaviartThomas(S, 2, variant="equispaced")',
'RaviartThomas(S, 3, variant="equispaced")',
"DiscontinuousRaviartThomas(T, 1)",
"DiscontinuousRaviartThomas(T, 2)",
"DiscontinuousRaviartThomas(T, 3)",
Expand All @@ -180,12 +180,12 @@ def __init__(self, a, b):
'BrezziDouglasMarini(S, 1, variant="integral(2)")',
'BrezziDouglasMarini(S, 2, variant="integral(3)")',
'BrezziDouglasMarini(S, 3, variant="integral(4)")',
'BrezziDouglasMarini(T, 1, variant="point")',
'BrezziDouglasMarini(T, 2, variant="point")',
'BrezziDouglasMarini(T, 3, variant="point")',
'BrezziDouglasMarini(S, 1, variant="point")',
'BrezziDouglasMarini(S, 2, variant="point")',
'BrezziDouglasMarini(S, 3, variant="point")',
'BrezziDouglasMarini(T, 1, variant="equispaced")',
'BrezziDouglasMarini(T, 2, variant="equispaced")',
'BrezziDouglasMarini(T, 3, variant="equispaced")',
'BrezziDouglasMarini(S, 1, variant="equispaced")',
'BrezziDouglasMarini(S, 2, variant="equispaced")',
'BrezziDouglasMarini(S, 3, variant="equispaced")',
"Nedelec(T, 1)",
"Nedelec(T, 2)",
"Nedelec(T, 3)",
Expand All @@ -204,12 +204,12 @@ def __init__(self, a, b):
'Nedelec(S, 1, variant="integral(2)")',
'Nedelec(S, 2, variant="integral(3)")',
'Nedelec(S, 3, variant="integral(4)")',
'Nedelec(T, 1, variant="point")',
'Nedelec(T, 2, variant="point")',
'Nedelec(T, 3, variant="point")',
'Nedelec(S, 1, variant="point")',
'Nedelec(S, 2, variant="point")',
'Nedelec(S, 3, variant="point")',
'Nedelec(T, 1, variant="equispaced")',
'Nedelec(T, 2, variant="equispaced")',
'Nedelec(T, 3, variant="equispaced")',
'Nedelec(S, 1, variant="equispaced")',
'Nedelec(S, 2, variant="equispaced")',
'Nedelec(S, 3, variant="equispaced")',
"NedelecSecondKind(T, 1)",
"NedelecSecondKind(T, 2)",
"NedelecSecondKind(T, 3)",
Expand All @@ -228,12 +228,12 @@ def __init__(self, a, b):
'NedelecSecondKind(S, 1, variant="integral(2)")',
'NedelecSecondKind(S, 2, variant="integral(3)")',
'NedelecSecondKind(S, 3, variant="integral(4)")',
'NedelecSecondKind(T, 1, variant="point")',
'NedelecSecondKind(T, 2, variant="point")',
'NedelecSecondKind(T, 3, variant="point")',
'NedelecSecondKind(S, 1, variant="point")',
'NedelecSecondKind(S, 2, variant="point")',
'NedelecSecondKind(S, 3, variant="point")',
'NedelecSecondKind(T, 1, variant="equispaced")',
'NedelecSecondKind(T, 2, variant="equispaced")',
'NedelecSecondKind(T, 3, variant="equispaced")',
'NedelecSecondKind(S, 1, variant="equispaced")',
'NedelecSecondKind(S, 2, variant="equispaced")',
'NedelecSecondKind(S, 3, variant="equispaced")',
"Regge(T, 0)",
"Regge(T, 1)",
"Regge(T, 2)",
Expand Down Expand Up @@ -546,7 +546,7 @@ def test_expansion_values(dim):
dpoints.append(tuple(2*np.array(alpha, dtype="d")/npoints-1))
rpoints.append(tuple(2*sympy.Rational(a, npoints)-1 for a in alpha))

n = 48
n = 20
eta = sympy.DeferredVector("eta")
Uvals = U.tabulate(n, dpoints)
if dim == 1:
Expand All @@ -561,8 +561,8 @@ def test_expansion_values(dim):
assert error < 1E-13
elif dim == 2:
idx = expansions.morton_index2
for p in range(n + 1):
q = n - p
for q in range(n + 1):
p = n - q
f = (sympy.jacobi_poly(p, 0, 0, eta[0]) *
sympy.jacobi_poly(q, 2*p+1, 0, eta[1]) * ((1 - eta[1])/2) ** p)
f *= sympy.sqrt((half + p) * (1 + p + q))
Expand Down

0 comments on commit 79a47d0

Please sign in to comment.