Hi Tuomas,
  If you are using hexahedra, a workaround is to build the horizontal part of RT1 using OuterProductElement and Hdiv. This has the advantage of doing the right thing on the sphere and over topography.

all the best
--Colin

On 20 May 2015 at 02:22, Tuomas Karna <tuomas.karna@gmail.com> wrote:
Hi all,

I'm implementing 3D momentum equation for horizontal velocity components (u,v) that live in RT1xP1 space. So the space is 3D but the field only has 2 components. I'm getting some errors from FFC related to missing quadrature rules:

File "/home/tuomas/.local/lib/python2.7/site-packages/FFC-1.6.0.dev0-py2.7-linux-x86_64.egg/ffc/quadrature/quadraturetransformerbase.py", line 1031, in _create_mapping_basis
    name, non_zeros, zeros, ones = self.name_map[name]
KeyError: 'FE1_C2'

If I interpret that correctly, FE1_C2 is a quad rule for the third component of the vector space, which doesn't exist in this case. If I replace the RT space with DG (dim=2) everything works. Is this a FFC issue or am I doing this wrong?

- Tuomas

---

from firedrake import *

mesh2d = UnitSquareMesh(10, 10)
layers = 6
mesh = ExtrudedMesh(mesh2d, layers, -1.0/layers)

#U = VectorFunctionSpace(mesh, 'DG', 1, vfamily='CG', vdegree=1, dim=2)
U = FunctionSpace(mesh, 'RT', 1, vfamily='CG', vdegree=1)
P1 = FunctionSpace(mesh, 'CG', 1, vfamily='CG', vdegree=1)

f = Function(P1)
f.interpolate(Expression('-x[0]'))
sol = Function(U)

test = TestFunction(U)

# sol + grad(f) = 0
F = inner(sol, test)*dx + (grad(f)[0]*test[0] + grad(f)[1]*test[1])*dx

solve(F == 0, sol)
print sol.dat.data.min(), sol.dat.data.max()


_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake