test_high_order_meshes.py

from firedrake import *
import numpy as np
import pytest


@pytest.mark.xfail(reason="Degree estimation wrong")
def test_high_order_coords():
    m = UnitTriangleMesh()
    V = VectorFunctionSpace(m, 'CG', 3)
    f = Function(V)
    f.project(m.coordinates)

    # Warp mesh so that the bottom triangle line is:
    # x(x - 1)(x + a) with a = 19/12.0
    def warp(x):
        return x * (x - 1)*(x + 19/12.0)

    f.dat.data[1, 1] = warp(1.0/3.0)
    f.dat.data[2, 1] = warp(2.0/3.0)

    mesh = Mesh(f)

    assert np.allclose(assemble(1*dx(domain=mesh)),
                       0.5 - (1.0/4.0 - (1 - 19.0/12.0)/3.0 - 19/24.0))