From 27589c03b44e3e28b3b74647520e69185cd07280 Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Fri, 28 Mar 2014 16:58:18 +0000 Subject: [PATCH 1/4] Fix output of VFS on interval --- firedrake/io.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/firedrake/io.py b/firedrake/io.py index 1a4471d588..980cacaaf4 100644 --- a/firedrake/io.py +++ b/firedrake/io.py @@ -277,8 +277,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) From 58a6e78ab9a1d2eaa792d266aa7017206d2982b5 Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Mon, 31 Mar 2014 10:50:19 +0100 Subject: [PATCH 2/4] Fix output of HDiv/HCurl fields The test for if they were CG or not was wrong. --- firedrake/io.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/firedrake/io.py b/firedrake/io.py index 980cacaaf4..091acec8f2 100644 --- a/firedrake/io.py +++ b/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': From ee63a1acbefb2840f6c2171be2105466c9d5bfe8 Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Fri, 28 Mar 2014 16:58:06 +0000 Subject: [PATCH 3/4] Fix DG coordinate field for periodic interval Ugly hack until we get fed it from DMPlex. Note that this is probably wrong when symbolic geometry stuff lands. --- firedrake/mesh.py | 29 ++++++++--------------------- 1 file changed, 8 insertions(+), 21 deletions(-) diff --git a/firedrake/mesh.py b/firedrake/mesh.py index 3f39357ad2..04d191fd63 100644 --- a/firedrake/mesh.py +++ b/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) From f7f7f905a1e27e2b753d2624c79ab6b95a5fb758 Mon Sep 17 00:00:00 2001 From: Lawrence Mitchell Date: Mon, 31 Mar 2014 11:54:50 +0100 Subject: [PATCH 4/4] Add some simple div-free tests for extruded intervals --- tests/extrusion/test_extrusion_interval.py | 75 +++++++++++++++++++ .../extrusion/test_extrusion_interval_area.py | 25 ------- 2 files changed, 75 insertions(+), 25 deletions(-) create mode 100644 tests/extrusion/test_extrusion_interval.py delete mode 100644 tests/extrusion/test_extrusion_interval_area.py diff --git a/tests/extrusion/test_extrusion_interval.py b/tests/extrusion/test_extrusion_interval.py new file mode 100644 index 0000000000..5c4c5d6193 --- /dev/null +++ b/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__)) diff --git a/tests/extrusion/test_extrusion_interval_area.py b/tests/extrusion/test_extrusion_interval_area.py deleted file mode 100644 index a43d02808e..0000000000 --- a/tests/extrusion/test_extrusion_interval_area.py +++ /dev/null @@ -1,25 +0,0 @@ -import pytest - -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 - -if __name__ == '__main__': - import os - pytest.main(os.path.abspath(__file__))