Re: [firedrake] 3d momentum equation with RT spaces
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()
Hi Andrew, Thanks for this, makes sense. So far I've used P1 in vertical because some operations were broken for discontinuous fields, but I'll give it a shot. - Tuomas On 05/20/2015 01:18 AM, Andrew McRae wrote:
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 <mailto: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
participants (2)
-
Andrew McRae
-
Tuomas Karna