enriched element bork

[dorugeber]

from firedrake import *

m = UnitTriangleMesh()
mesh = ExtrudedMesh(m, layers=3, layer_height=1.0)
V0 = FunctionSpace(mesh, "CG", 1, vfamily="CG", vdegree=1)

V1_a_horiz = FiniteElement("RT", "triangle", 1)
V1_a_vert = FiniteElement("CG", "interval", 1)
V1_a = HCurl(OuterProductElement(V1_a_horiz, V1_a_vert))

V1_b_horiz = FiniteElement("CG", "triangle", 2)
V1_b_vert = FiniteElement("DG", "interval", 0)
V1_b = HCurl(OuterProductElement(V1_b_horiz, V1_b_vert))

V1_elt = EnrichedElement(V1_a, V1_b)
V1 = FunctionSpace(mesh, V1_elt)

u = Function(V0)
u.interpolate(Expression("x[2]"))

v = TrialFunction(V1)
w = TestFunction(V1)
a = dot(v, w)*dx
L = dot(grad(u), w)*dx

v = Function(V1)
parms = {'snes_type': 'ksponly', 'ksp_type': 'preonly', 'pc_type': 'lu'}
solve(a == L, v, solver_parameters=parms)
print np.sort(v.dat.data) # expect 12 1s and 9 zeros

[dorugeber]

RHS might be correct [probably not in general], but LHS matrix is incorrect. Can show using eigenvalue decomposition, say.

[dorugeber]

Possible numbering issue.

First 6 dofs are from the RT part (horizontal), last 6 are from the P2 x P0 part (vertical)

In [13]: V1.cell_node_map().values
Out[13]: array([[14, 15,  7,  8,  2,  3,  0, 12, 19, 16,  9,  4]], dtype=int32)

In [14]: V1.offset + V1.cell_node_map().values
Out[14]: array([[16, 17,  9, 10,  4,  5,  1, 13, 20, 18, 11,  6]], dtype=int32)

Note that node 4 turns from a vertical to a horizontal dof (WRONG)

Looks like the numbering up the column goes "2, 4, 3, 6, 5", so although the offsets are correct, they don't work in the first layer, since node 2 + offset 2 lands at 4 instead of 3.

[wence-]

We think this is a FIAT numbering bug:

P2 = FunctionSpace(mesh, 'CG', 2, vfamily='CG', vdegree=2)

P2.fiat_element.entity_dofs()
{(0, 0): {0: [0], 1: [1], 2: [3], 3: [4], 4: [6], 5: [7]},
 (0, 1): {0: [2], 1: [5], 2: [8]},
 (1, 0): {0: [9], 1: [10], 2: [12], 3: [13], 4: [15], 5: [16]},
 (1, 1): {0: [11], 1: [14], 2: [17]},
 (2, 0): {0: [], 1: []},
 (2, 1): {0: []}}

P2.fiat_element.flattened_element().entity_dofs()

{0: {0: [0, 2, 1], 1: [3, 5, 4], 2: [6, 8, 7]},
 1: {0: [9, 11, 10], 1: [12, 14, 13], 2: [15, 17, 16]},
 2: {0: []}}

V1.fiat_element.entity_dofs()

{(0, 0): {0: [], 1: [], 2: [], 3: [], 4: [], 5: []},
 (0, 1): {0: [6], 1: [7], 2: [8]},
 (1, 0): {0: [0], 1: [1], 2: [2], 3: [3], 4: [4], 5: [5]},
 (1, 1): {0: [9], 1: [10], 2: [11]},
 (2, 0): {0: [], 1: []},
 (2, 1): {0: []}}

V1.fiat_element.flattened_element().entity_dofs()
{0: {0: [6], 1: [7], 2: [8]},
 1: {0: [0, 1, 9], 1: [2, 3, 10], 2: [4, 5, 11]},
 2: {0: []}}

Note how in the P2 case the dofs from (1,1) have been interleaved into the dofs from (1, 0), whereas in the enriched case they have not. In global numbering we hand things out in this order, which is why things are wrong.

In particular, the entity_dofs of the V1 flattened element should be (we think).

{0: {0: [6], 1: [7], 2: [8]},
 1: {0: [0, 9, 1], 1: [2, 10, 3], 2: [4, 11, 5]},
 2: {0: []}}

then everything would fall out.

[dorugeber]

Oh god... yes, TFE is currently hacked to output dofs from bottom to top, as Doru requested originally. I'll think about how this can be done for enriched elements.

[dorugeber]

What is the minimum that needs to be done? Full sorting is probably overkill, and unnecessarily hard (eg sorting the internal dofs of P3 and P4). I will have a think.

[wence-]

I think we just need shared ones to come out numbered after non-shared, full bottom to top sorting would be nice. Is it just a case of interchanging the order in which you concatenate the entities of a given dimension?

[dorugeber]

Yeah, at the moment I think I generate the flattened dofs by concatenating the flattened dofs of the constituent elements. I think I know how to unhack and fix this.

[dorugeber]

(it involves unhacking the actual flattening process)

[dorugeber]

Mesdames et messieurs, je présente un Helmholtz mixte (avec un peu de superconvergence)

"""Tests for mixed Helmholtz convergence on extruded meshes"""
import numpy as np

from firedrake import *

for ii in [1, 2, 3, 4, 5]:
    m = UnitSquareMesh(2**ii, 2**ii)
    mesh = ExtrudedMesh(m, layers=(2**ii + 1), layer_height = 1.0/2**ii)

    V2_a_horiz = FiniteElement("BDM", "triangle", 1)
    V2_a_vert = FiniteElement("DG", "interval", 0)
    V2_a = HDiv(OuterProductElement(V2_a_horiz, V2_a_vert))

    V2_b_horiz = FiniteElement("DG", "triangle", 0)
    V2_b_vert = FiniteElement("CG", "interval", 1)
    V2_b = HDiv(OuterProductElement(V2_b_horiz, V2_b_vert))

    V2_elt = EnrichedElement(V2_a, V2_b)
    V2 = FunctionSpace(mesh, V2_elt)

    V3 = FunctionSpace(mesh, "DG", 0, vfamily="DG", vdegree=0)

    f = Function(V3)
    exact = Function(V3)
    f.interpolate(Expression("(1+12*pi*pi)*sin(x[0]*pi*2)*sin(x[1]*pi*2)*sin(x[2]*pi*2)"))
    exact.interpolate(Expression("sin(x[0]*pi*2)*sin(x[1]*pi*2)*sin(x[2]*pi*2)"))

    W = V2 * V3
    u, p = TrialFunctions(W)
    v, q = TestFunctions(W)
    a = (p*q - q*div(u) + dot(v, u) + div(v)*p)*dx
    L = f*q*dx

    out = Function(W)
    solve(a == L, out, solver_parameters={'pc_type': 'fieldsplit',
                                          'pc_fieldsplit_type': 'schur',
                                          'ksp_type': 'cg',
                                          'pc_fieldsplit_schur_fact_type': 'FULL',
                                          'fieldsplit_0_ksp_type': 'cg',
                                          'fieldsplit_1_ksp_type': 'cg'})
    print "L_inf norm: " + str(max(abs(out.dat[1].data - exact.dat.data)))
    print "L2 norm: " + str(sqrt(assemble((out[3]-exact)*(out[3]-exact)*dx)))

[doru1004]

Ca c'est magnifique (magnifeec!?) !!!

[kynan]

In which revision was this fixed?

[wence-]

It was a FIAT bug, fixed in FIAT#9 (enriched) merge.