Hi Tuomas,

Are you sure you want to use RT1xP1, as opposed to RT1xP0[DG]? This would be the usual element for 'mimetic' discretisations.

The code for RT1xP0 would be [possibly with typos...]

---------------------------
from firedrake import *

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

Uh_elt = FiniteElement('RT', triangle, 1)
Uv_elt = FiniteElement('DG', interval, 0)
U_elt = HDiv(OuterProductElement(Uh_elt, Uv_elt))
U = FunctionSpace(mesh, U_elt)

P1 = FunctionSpace(mesh, 'CG', 1) # automatically gives the appropriate P1 x P1 element

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

test = TestFunction(U)

# sol + grad(f) = 0
F = inner(sol, test)*dx + dot(grad(f), test)*dx

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

FIAT will throw an error if you try to do this with RT1xP1, however (though it will let you build a curl-conforming 'edge' element by using HCurl). If you genuinely need RT1xP1 then I can probably turn off some FIAT validation so that it can be used. I think it's more likely that you want RT1xP0, however, and a corresponding vertical velocity element of P0xP1, HDiv'd.

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()