Projection into DG1 x CG1 space looks wrong

[wence-]

The following code:

from firedrake import *

L = 2.0
H = 1.0
nlayers = 100
dz = H/nlayers

m = PeriodicIntervalMesh(100, length=L)
mesh = ExtrudedMesh(m, layers=nlayers, layer_height=dz)

#1D spaces in horizontal
IUh0 = FiniteElement("CG", "interval", 2)
IUh1 = FiniteElement("DG", "interval", 1)

#1D spaces in vertical
IUv0 = FiniteElement("CG", "interval", 1)
IUv1 = FiniteElement("DG", "interval", 0)

#Temperature space
V1v_elt = OuterProductElement(IUh1, IUv0)
V1v = FunctionSpace(mesh, V1v_elt)

#Solution variables
oldT = project(Expression(("exp(-((x[0]-1.0)*(x[0]-1.0) + (x[1]-0.5)*(x[1]-0.5))/0.04)")),V1v)

file = File("temp.pvd")

file << oldT

produces a field where the slopes in each element go the wrong way. This can most easily be observed by opening the file in Paraview and selecting Filters -> Warp by Scalar, then clicking Apply.

Is something wrong with the numbering perhaps?

[dorugeber]

Can you wrap your code in ``` marks? Otherwise, some characters (notably '*') are interpreted as formatting, so your code comes out with random italics, and lacking characters.
Edit: done it.

[dorugeber]

Probably a dmplex thing, one sec....

[dorugeber]

Nope, my first idea didn't work.

[dorugeber]

Note that dmplex changed the interval numbering from
0----1 2----3 4----5
to
1----0 3----2 5----4
which must be related...

[colinjcotter]

I'm going to check this for DG as well, and will make a failing test.

[colinjcotter]

Some further experimentation indicates that this is probably a bug with the vtu output, not FFC.

from firedrake import *

L = 1.0
H = 1.0
nlayers = 5
dz = H/nlayers

m = PeriodicIntervalMesh(5, length=1.0)
mesh = ExtrudedMesh(m, layers=nlayers, layer_height=1.0/5.0)

#1D spaces in horizontal
IUh0 = FiniteElement("CG", "interval", 2)
IUh1 = FiniteElement("DG", "interval", 1)

#1D spaces in vertical
IUv0 = FiniteElement("CG", "interval", 1)
IUv1 = FiniteElement("DG", "interval", 0)

#P0 space
V0_elt = OuterProductElement(IUh0, IUv0)
V0 = FunctionSpace(mesh, V0_elt)

#DG space
V2_elt = OuterProductElement(IUh1, IUv1)
V2 = FunctionSpace(mesh, V2_elt)

#Solution variables
T = project(Expression(("x[0]")),V0)

file = File("temp.pvd")

file << T

[wence-]

There are three bugs here.

1. the Project function does the wrong thing. It should interpolate Expressions into L2 not H1 before projecting
2. The coordinate field on periodic meshes is incorrect given the new global numbering
3. Numbering on periodic intervals is somehow broken (at least for interpolation):

from firedrake import *


m = PeriodicUnitIntervalMesh(4)

V = FunctionSpace(m, 'CG', 1)

print m.coordinates.dat.data

u = Function(V)


# Interpolate does this:
to_element = V.fiat_element

to_pts = []
for dual in to_element.dual_basis():
    to_pts.append(dual.pt_dict.keys()[0])

coords = m.coordinates
coords_space = coords.function_space()
coords_element = coords_space.fiat_element
X = coords_element.tabulate(0, to_pts).values()[0]

for cell in range(4):
    for k in range(2):
        print cell, u.cell_node_map().values[cell][k], sum([X.T[k][i] * coords.dat.data[coords.cell_node_map().values[cell]][i] for i in range(2)])
        u.dat.data[u.cell_node_map().values[cell][k]] = sum([X.T[k][i] * coords.dat.data[coords.cell_node_map().values[cell]][i] for i in range(2)])

print u.dat.data

I expect to see [1, 0.25, 0.5, 0.75] at the end

[wence-]

Actual output:

[ 0.    0.25  0.25  0.5   0.5   0.75  0.75  1.  ]
0 1 0.0
0 0 0.25
1 2 0.25
1 1 0.5
2 3 0.5
2 2 0.75
3 3 0.75
3 0 1.0
[ 1.    0.5   0.75  0.75]

This is after I have fixed bug 2.

[wence-]

Michael, this seems to be a hangover from the move to dmplex due to facet numbering. I'm lost in a world of cones, all alike. Can you take a look?

[mlange05]

OK, the cell numbering in the UnitInterval Plex does not agree with the cell numbering of the periodic coordinates:

Plex does 0, 1, ..., c, but the periodic coordinates are written for 0, c, 2, 3, ..., c-1.

[wence-]

I think I have a fix for this. Will push when I get home. It's really ugly though. Colin, once I update, can you have another go with various pieces , or suggest done better ways if testing for breakages other than computing areas. I can think of some simple minded things computing interior facet integrals dotted against facet normals but there must be others.

[wence-]

fix/periodic does this, I think. But I need to think harder, I'm not really sure what's going on.

[wence-]

I /think/ this is fixed by this branch. But I would like some very careful eyes over this bit of the code. Here's a summary of various bits and pieces I have noticed.

1. dmplex doesn't distinguish between facets and vertices in 1d, and does some magic vertex renumbering, I think to get facet normal tests to pass. I /think/ this should be undone, and replaced with appropriate construction of "reversed" facet maps for 1d and extruded 1d.
2. As a corollary to this, the DG coordinate field needs all its values feeding in "backwards". And the last element is "reversed". I'm sort of not happy with this hack.
3. The tests in this area are pretty weak. I'm not yet convinced that interior facets are correct for periodic intervals. I have added two extruded tests (on periodic and non-periodic) checking that a field I claim is divergence-free actually is. Please ensure these are doing the right thing.

[dham]

I believe the hack to be correct, if aesthetically vile.

[dorugeber]

Is there a good reason why the last cell is reversed?

[dorugeber]

And are you implying some new/modified FFC code snippets for the facet normals? :)

[wence-]

All existing facet normal code still passes with this change. So one might need some new code, but no existing code needs changing. Note that there are basically no tests of facets on periodic, so things may still be broken there.

The last cell needs to be reversed so that the coordinate field lines up with the numbering. That's not a good reason, I can draw some pictures tomorrow if you'd like.

[dorugeber]

Believeable.

[wence-]

Merged and hopefully fixed. Please, however, keep a keen eye on the periodic pieces, especially facets when they arrive, because these are still only very lightly tested.