Skip to content

Commit

Permalink
Revert "Merged in branch 'toriesz' (pull request #54)"
Browse files Browse the repository at this point in the history 
This reverts commit 1375329, reversing
changes made to daacf42.

This breaks the Bell element, and also FacetBubble.
  • Loading branch information
@wence-
wence- committed Feb 19, 2019
1 parent 1375329 commit 54c4a4a
Show file tree
Hide file tree
Showing 7 changed files with 102 additions and 189 deletions.
73 changes: 4 additions & 69 deletions FIAT/dual_set.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,9 +16,6 @@
# along with FIAT. If not, see <http://www.gnu.org/licenses/>.

import numpy
import collections

from FIAT import polynomial_set


class DualSet(object):
Expand Down Expand Up @@ -53,79 +50,17 @@ def get_reference_element(self):
return self.ref_el

def to_riesz(self, poly_set):
"""This method gives the action of the entire dual set
on each member of the expansion set underlying poly_set.
Then, applying the linear functionals of the dual set to an
arbitrary polynomial in poly_set is accomplished by matrix
multiplication."""

# This rather technical code queries the low-level information
# for each functional to find out where it evaluates its
# inputs and/or their derivatives. Then, it tabulates the
# expansion set one time for all the function values and
# another for all of the derivatives. This circumvents
# needing to call the to_riesz method of each functional and
# also limits the number of different calls to tabulate.

tshape = self.nodes[0].target_shape
num_nodes = len(self.nodes)
es = poly_set.get_expansion_set()
ed = poly_set.get_embedded_degree()
num_exp = es.get_num_members(poly_set.get_embedded_degree())

riesz_shape = tuple([num_nodes] + list(tshape) + [num_exp])

result = numpy.zeros(riesz_shape, "d")

# let's amalgamate the pt_dict and deriv_dicts of all the
# functionals so we can tabulate the basis functions twice only
# (once on pts and once on derivatives)

# Need: dictionary mapping pts to which functionals they come from
pts_to_ells = collections.OrderedDict()
dpts_to_ells = collections.OrderedDict()

for i, ell in enumerate(self.nodes):
for pt in ell.pt_dict:
pts_to_ells.setdefault(pt, []).append(i)

for i, ell in enumerate(self.nodes):
for pt in ell.deriv_dict:
dpts_to_ells.setdefault(pt, []).append(i)

# Now tabulate
pts = list(pts_to_ells.keys())
expansion_values = es.tabulate(ed, pts)

for j, pt in enumerate(pts):
which_ells = pts_to_ells[pt]

for k in which_ells:
pt_dict = self.nodes[k].pt_dict
wc_list = pt_dict[pt]

for i in range(num_exp):
for (w, c) in wc_list:
result[k][c][i] += w*expansion_values[i, j]

max_deriv_order = max([ell.max_deriv_order for ell in self.nodes])
if max_deriv_order > 0:
dpts = list(dpts_to_ells.keys())
# It's easiest/most efficient to get derivatives of the
# expansion set through the polynomial set interface.
# This is creating a short-lived set to do just this.
expansion = polynomial_set.ONPolynomialSet(self.ref_el, ed)
dexpansion_values = expansion.tabulate(dpts, max_deriv_order)

for j, pt in enumerate(dpts):
which_ells = dpts_to_ells[pt]

for k in which_ells:
dpt_dict = self.nodes[k].deriv_dict
wac_list = dpt_dict[pt]
self.mat = numpy.zeros(riesz_shape, "d")

for i in range(num_exp):
for (w, alpha, c) in wac_list:
result[k][c][i] += w*dexpansion_values[alpha][i, j]
for i in range(len(self.nodes)):
self.mat[i][:] = self.nodes[i].to_riesz(poly_set)

return result
return self.mat
90 changes: 43 additions & 47 deletions FIAT/functional.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,6 @@
import numpy
import sympy

from FIAT import polynomial_set


def index_iterator(shp):
"""Constructs a generator iterating over all indices in
Expand All @@ -44,6 +42,11 @@ def index_iterator(shp):
for foo in index_iterator(shp_foo):
yield [i] + foo

# also put in a "jet_dict" that maps
# pt --> {wt, multiindex, comp}
# the multiindex is an iterable of nonnegative
# integers


class Functional(object):
"""Class implementing an abstract functional.
Expand Down Expand Up @@ -91,6 +94,7 @@ def get_type_tag(self):
normal component, which is probably handy for clients of FIAT"""
return self.functional_type

# overload me in subclasses to make life easier!!
def to_riesz(self, poly_set):
"""Constructs an array representation of the functional over
the base of the given polynomial_set so that f(phi) for any
Expand All @@ -99,20 +103,17 @@ def to_riesz(self, poly_set):
ed = poly_set.get_embedded_degree()
pt_dict = self.get_point_dict()

nbf = poly_set.get_num_members()

pts = list(pt_dict.keys())

# bfs is matrix that is pdim rows by num_pts cols
# where pdim is the polynomial dimension

bfs = es.tabulate(ed, pts)
npts = len(pts)

result = numpy.zeros(poly_set.coeffs.shape[1:], "d")

# loop over points
for j in range(npts):
for j in range(len(pts)):
pt_cur = pts[j]
wc_list = pt_dict[pt_cur]

Expand All @@ -122,23 +123,7 @@ def to_riesz(self, poly_set):
result[c][i] += w * bfs[i, j]

if self.deriv_dict:
dpt_dict = self.deriv_dict

# this makes things quicker since it uses dmats after
# instantiation
es_foo = polynomial_set.ONPolynomialSet(self.ref_el, ed)
dpts = list(dpt_dict.keys())

# figure out the maximum order of derivative being taken
dbfs = es_foo.tabulate(dpts, self.max_deriv_order)

ndpts = len(dpts)
for j in range(ndpts):
dpt_cur = dpts[j]
wac_list = dpt_dict[dpt_cur]
for i in range(nbf):
for (w, alpha, c) in wac_list:
result[c][i] += w * dbfs[alpha][i, j]
raise NotImplementedError("Generic to_riesz implementation does not support derivatives")

return result

Expand Down Expand Up @@ -187,8 +172,8 @@ class PointDerivative(Functional):
functions at a particular point x."""

def __init__(self, ref_el, x, alpha):
dpt_dict = {x: [(1.0, tuple(alpha), tuple())]}
self.alpha = tuple(alpha)
dpt_dict = {x: [(1.0, alpha, tuple())]}
self.alpha = alpha
self.order = sum(self.alpha)

Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointDeriv")
Expand All @@ -206,6 +191,25 @@ def __call__(self, fn):

return sympy.diff(fn(X), *dvars).evalf(subs=dict(zip(dX, x)))

def to_riesz(self, poly_set):
x = list(self.deriv_dict.keys())[0]

X = sympy.DeferredVector('x')
dx = numpy.asarray([X[i] for i in range(len(x))])

es = poly_set.get_expansion_set()
ed = poly_set.get_embedded_degree()

bfs = es.tabulate(ed, [dx])[:, 0]

# Expand the multi-index as a series of variables to
# differentiate with respect to.
dvars = tuple(d for d, a in zip(dx, self.alpha)
for count in range(a))

return numpy.asarray([sympy.lambdify(X, sympy.diff(b, *dvars))(x)
for b in bfs])


class PointNormalDerivative(Functional):

Expand All @@ -219,35 +223,27 @@ def __init__(self, ref_el, facet_no, pt):
alpha = [0] * sd
alpha[i] = 1
alphas.append(alpha)
dpt_dict = {pt: [(n[i], tuple(alphas[i]), tuple()) for i in range(sd)]}
dpt_dict = {pt: [(n[i], alphas[i], tuple()) for i in range(sd)]}

Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv")

def to_riesz(self, poly_set):
x = list(self.deriv_dict.keys())[0]

class PointNormalSecondDerivative(Functional):
X = sympy.DeferredVector('x')
dx = numpy.asarray([X[i] for i in range(len(x))])

def __init__(self, ref_el, facet_no, pt):
n = ref_el.compute_normal(facet_no)
self.n = n
sd = ref_el.get_spatial_dimension()
tau = numpy.zeros((sd*(sd+1)//2,))
es = poly_set.get_expansion_set()
ed = poly_set.get_embedded_degree()

alphas = []
cur = 0
for i in range(sd):
for j in range(i, sd):
alpha = [0] * sd
alpha[i] += 1
alpha[j] += 1
alphas.append(tuple(alpha))
tau[cur] = n[i]*n[j]
cur += 1

self.tau = tau
self.alphas = alphas
dpt_dict = {pt: [(n[i], alphas[i], tuple()) for i in range(sd)]}
bfs = es.tabulate(ed, [dx])[:, 0]

Functional.__init__(self, ref_el, tuple(), {}, dpt_dict, "PointNormalDeriv")
# We need the gradient dotted with the normal.
return numpy.asarray(
[sympy.lambdify(
X, sum([sympy.diff(b, dxi)*ni
for dxi, ni in zip(dx, self.n)]))(x)
for b in bfs])


class IntegralMoment(Functional):
Expand Down
4 changes: 2 additions & 2 deletions bitbucket-pipelines.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,9 +6,9 @@ pipelines:
script:
- apt-get update
- apt-get install -y
git python3-minimal python3-pip python3-setuptools python3-wheel
git python3-minimal python3-pip python3-pytest python3-setuptools python3-wheel
--no-install-recommends
- pip3 install flake8 pytest --upgrade
- pip3 install flake8
- pip3 install .
- python3 -m flake8
- DATA_REPO_GIT="" python3 -m pytest -v test/
2 changes: 1 addition & 1 deletion test/regression/test_regression.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -277,7 +277,7 @@ def _perform_test(family, dim, degree, reference_table):
reference = json.load(open(filename, "r"), object_hook=json_numpy_obj_hook)
for test_case in test_cases:
family, dim, degree = test_case
_perform_test(family, dim, degree, reference[str(test_case)])
yield _perform_test, family, dim, degree, reference[str(test_case)]

# Update references if missing
except (IOError, KeyError) as e:
Expand Down
64 changes: 32 additions & 32 deletions test/unit/test_fiat.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -203,42 +203,42 @@ def __init__(self, a, b):

# MixedElement made of nodal elements should be nodal, but its API
# is currently just broken.
pytest.param("MixedElement(["
" DiscontinuousLagrange(T, 1),"
" RaviartThomas(T, 2)"
"])", marks=xfail_impl),
xfail_impl("MixedElement(["
" DiscontinuousLagrange(T, 1),"
" RaviartThomas(T, 2)"
"])"),

# Following element do not bother implementing get_nodal_basis
# so the test would need to be rewritten using tabulate
pytest.param("TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2))", marks=xfail_impl),
pytest.param("Hdiv(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))", marks=xfail_impl),
pytest.param("Hcurl(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))", marks=xfail_impl),
pytest.param("HDivTrace(T, 1)", marks=xfail_impl),
pytest.param("EnrichedElement("
"Hdiv(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), "
"Hdiv(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
")", marks=xfail_impl),
pytest.param("EnrichedElement("
"Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), "
"Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
")", marks=xfail_impl),
xfail_impl("TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2))"),
xfail_impl("Hdiv(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"),
xfail_impl("Hcurl(TensorProductElement(DiscontinuousLagrange(I, 1), Lagrange(I, 2)))"),
xfail_impl("HDivTrace(T, 1)"),
xfail_impl("EnrichedElement("
"Hdiv(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), "
"Hdiv(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
")"),
xfail_impl("EnrichedElement("
"Hcurl(TensorProductElement(Lagrange(I, 1), DiscontinuousLagrange(I, 0))), "
"Hcurl(TensorProductElement(DiscontinuousLagrange(I, 0), Lagrange(I, 1)))"
")"),
# Following elements are checked using tabulate
pytest.param("HDivTrace(T, 0)", marks=xfail_impl),
pytest.param("HDivTrace(T, 1)", marks=xfail_impl),
pytest.param("HDivTrace(T, 2)", marks=xfail_impl),
pytest.param("HDivTrace(T, 3)", marks=xfail_impl),
pytest.param("HDivTrace(S, 0)", marks=xfail_impl),
pytest.param("HDivTrace(S, 1)", marks=xfail_impl),
pytest.param("HDivTrace(S, 2)", marks=xfail_impl),
pytest.param("HDivTrace(S, 3)", marks=xfail_impl),
pytest.param("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))", marks=xfail_impl),
pytest.param("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))", marks=xfail_impl),
pytest.param("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))", marks=xfail_impl),
pytest.param("TensorProductElement(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)), Lagrange(I, 2))", marks=xfail_impl),
pytest.param("FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)))", marks=xfail_impl),
pytest.param("FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)))", marks=xfail_impl),
pytest.param("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))", marks=xfail_impl),
pytest.param("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))", marks=xfail_impl),
xfail_impl("HDivTrace(T, 0)"),
xfail_impl("HDivTrace(T, 1)"),
xfail_impl("HDivTrace(T, 2)"),
xfail_impl("HDivTrace(T, 3)"),
xfail_impl("HDivTrace(S, 0)"),
xfail_impl("HDivTrace(S, 1)"),
xfail_impl("HDivTrace(S, 2)"),
xfail_impl("HDivTrace(S, 3)"),
xfail_impl("TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))"),
xfail_impl("TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))"),
xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)), Lagrange(I, 1))"),
xfail_impl("TensorProductElement(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)), Lagrange(I, 2))"),
xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1)))"),
xfail_impl("FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2)))"),
xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 1), Lagrange(I, 1))), Lagrange(I, 1)))"),
xfail_impl("FlattenedDimensions(TensorProductElement(FlattenedDimensions(TensorProductElement(Lagrange(I, 2), Lagrange(I, 2))), Lagrange(I, 2)))"),
]


Expand Down
Loading

0 comments on commit 54c4a4a

Please sign in to comment.