Skip to content

Commit

Permalink
Merge pull request #205 from firedrakeproject/wence/fix/cellsize-facets
Browse files Browse the repository at this point in the history
Fix indexing of cell_size on interior face terms
  • Loading branch information
wence- authored Feb 5, 2020
2 parents 9bd9ec7 + e0c498b commit 6a0c1b7
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 3 deletions.
33 changes: 33 additions & 0 deletions tests/test_tsfc_204.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,33 @@
from tsfc import compile_form
from ufl import (BrokenElement, Coefficient, FacetNormal, FiniteElement,
FunctionSpace, Mesh, MixedElement, VectorElement, as_matrix,
dot, dS, ds, dx, facet, grad, inner, outer, split, triangle)


def test_physically_mapped_facet():
mesh = Mesh(VectorElement("P", triangle, 1))

# set up variational problem
U = FiniteElement("Morley", mesh.ufl_cell(), 2)
V = FiniteElement("P", mesh.ufl_cell(), 1)
R = FiniteElement("P", mesh.ufl_cell(), 1)
Vv = VectorElement(BrokenElement(V))
Qhat = VectorElement(BrokenElement(V[facet]))
Vhat = VectorElement(V[facet])
Z = FunctionSpace(mesh, MixedElement(U, Vv, Qhat, Vhat, R))

z = Coefficient(Z)
u, d, qhat, dhat, lam = split(z)

s = FacetNormal(mesh)
trans = as_matrix([[1, 0], [0, 1]])
mat = trans*grad(grad(u))*trans + outer(d, d) * u
J = (u**2*dx +
u**3*dx +
u**4*dx +
inner(mat, mat)*dx +
inner(grad(d), grad(d))*dx +
dot(s, d)**2*ds)
L_match = inner(qhat, dhat - d)
L = J + inner(lam, inner(d, d)-1)*dx + (L_match('+') + L_match('-'))*dS + L_match*ds
compile_form(L)
7 changes: 4 additions & 3 deletions tsfc/kernel_interface/common.py
Original file line number Diff line number Diff line change
Expand Up @@ -62,9 +62,10 @@ def cell_orientation(self, restriction):
def cell_size(self, restriction):
if not hasattr(self, "_cell_sizes"):
raise RuntimeError("Haven't called set_cell_sizes")
f = {None: (), '+': (0, ), '-': (1, )}[restriction]
# cell_sizes expression must have been set up by now.
return gem.partial_indexed(self._cell_sizes, f)
if self.interior_facet:
return self._cell_sizes[{'+': 0, '-': 1}[restriction]]
else:
return self._cell_sizes

def entity_number(self, restriction):
"""Facet or vertex number as a GEM index."""
Expand Down

0 comments on commit 6a0c1b7

Please sign in to comment.