Projection into DG1 x CG1 space looks wrong

firedrake/io.py

@@ -137,6 +137,9 @@ def __lshift__(self, data):
             function = data
 
         def is_family1(e, family):
+            import ufl.finiteelement.hdivcurl as hc
+            if isinstance(e, (hc.HDiv, hc.HCurl)):
+                return False
             if e.family() == 'OuterProductElement':
                 if e.degree() == (1, 1):
                     if e._A.family() == family \
@@ -147,6 +150,9 @@ def is_family1(e, family):
             return False
 
         def is_cgN(e):
+            import ufl.finiteelement.hdivcurl as hc
+            if isinstance(e, (hc.HDiv, hc.HCurl)):
+                return False
             if e.family() == 'OuterProductElement':
                 if e._A.family() == 'Lagrange' \
                    and e._B.family() == 'Lagrange':
@@ -277,8 +283,11 @@ def is_cgN(e):
         if isinstance(output.function_space(), VectorFunctionSpace):
             tmp = output.dat.data_ro_with_halos
             vdata = [None]*3
-            for i in range(output.dat.dim[0]):
-                vdata[i] = tmp[:, i].flatten()
+            if output.dat.dim[0] == 1:
+                vdata[0] = tmp.flatten()
+            else:
+                for i in range(output.dat.dim[0]):
+                    vdata[i] = tmp[:, i].flatten()
             for i in range(output.dat.dim[0], 3):
                 vdata[i] = np.zeros_like(vdata[0])
             data = tuple(vdata)

firedrake/mesh.py

@@ -302,28 +302,15 @@ def __init__(self, ncells, length):
         coordvec = PETSc.Vec().createWithArray(coords, size=size)
         dmplex.setCoordinatesLocal(coordvec)
 
-        # Coordinate values need to be replaced by the appropriate
-        # DG coordinate field.
         dx = length / ncells
-        # Two per cell
-        coords = np.empty(2 * ncells, dtype=float)
-        # For an interval
-        #
-        # 0---1---2---3 ... n-1---n
-        # |                       |
-        # `-----------------------'
-        #
-        # The element (0,1) is numbered first
-        coords[0] = 0.0
-        coords[1] = dx
-        # Then the element (n, 0)
-        coords[2] = length
-        coords[3] = length - dx
-        # Then the rest in order (1, 2), (2, 3) ... (n-1, n)
-        if len(coords) > 4:
-            coords[4] = dx
-            coords[5:] = np.repeat(np.arange(dx * 2, length - dx + dx*0.01, dx), 2)[:-1]
-
+        # HACK ALERT!
+        # Almost certainly not right when symbolic geometry stuff lands.
+        # Hopefully DMPlex will eventually give us a DG coordinate
+        # field.  Until then, we build one by hand.
+        coords = np.dstack((np.arange(dx, length + dx*0.01, dx),
+                            np.arange(0, length - dx*0.01, dx))).flatten()
+        # Last cell is back to front.
+        coords[-2:] = coords[-2:][::-1]
         Mesh.__init__(self, self.name, plex=dmplex,
                       periodic_coords=coords)
 

tests/extrusion/test_extrusion_interval.py

@@ -0,0 +1,75 @@
+import pytest
+import numpy as np
+from firedrake import *
+
+
+def integrate_one(intervals):
+    m = UnitIntervalMesh(intervals)
+    layers = intervals
+    mesh = ExtrudedMesh(m, layers, layer_height=1.0 / layers)
+
+    V = FunctionSpace(mesh, 'CG', 1)
+
+    u = Function(V)
+
+    u.interpolate(Expression("1"))
+
+    return assemble(u * dx)
+
+
+def test_unit_interval():
+    assert abs(integrate_one(5) - 1) < 1e-12
+
+
+def test_interval_div_free():
+    m = UnitIntervalMesh(50)
+    mesh = ExtrudedMesh(m, 50)
+
+    V = VectorFunctionSpace(mesh, 'CG', 3)
+
+    u = Function(V)
+
+    u.interpolate(Expression(('x[0]*x[0]*x[1]', '-x[0]*x[1]*x[1]')))
+
+    # u is pointwise divergence free, so the integral should also be
+    # div-free.
+    assert np.allclose(assemble(div(u)*dx), 0)
+
+    L2 = FunctionSpace(mesh, 'DG', 2)
+
+    v = TestFunction(L2)
+
+    f = assemble(div(u)*v*dx)
+
+    # Check pointwise div-free
+    assert np.allclose(f.dat.data, 0)
+
+
+def test_periodic_interval_div_free():
+    m = PeriodicUnitIntervalMesh(50)
+    mesh = ExtrudedMesh(m, 50)
+
+    V = VectorFunctionSpace(mesh, 'CG', 3)
+
+    u = Function(V)
+
+    u.interpolate(Expression(('sin(2*pi*x[0])',
+                              '-2*pi*x[1]*cos(2*pi*x[0])')))
+
+    # u is pointwise divergence free, so the integral should also be
+    # div-free.
+    assert np.allclose(assemble(div(u)*dx), 0)
+
+    L2 = FunctionSpace(mesh, 'DG', 2)
+
+    v = TestFunction(L2)
+
+    f = assemble(div(u)*v*dx)
+
+    # Check pointwise div-free
+    assert np.allclose(f.dat.data, 0)
+
+
+if __name__ == '__main__':
+    import os
+    pytest.main(os.path.abspath(__file__))