ufc::finite_element::evaluate_reference_basis_derivatives returns wrong values

[Prof Garth Wells]

The FFC generated UFC function ufc::finite_element::evaluate_reference_basis_derivatives is returning wrong results. For a P1 basis, I would expect all derivatives to be +1/-1, or 0. However, in 2D the function returns some two, and in 3D some threes.

Code to reproduce attached.

[Prof Garth Wells]

Change 6ca3182f8b76a3925a5db67d433732fec9ff3df5 might be related.

[Jan Blechta]

I would expect all derivatives to be +1/-1

I would expect also 0.5**-0.5 on reference triangle.

[Chris Richardson]

Not related to that change. The code for evaluate_reference_basis_derivatives is in ffc/uflacs/backends/ufc/evalderivs.py

[Prof Garth Wells]

@blechta Why would you expect 0.5**-0.5? It's P1 on the triangle reference triangle (0, 0), (1, 0), (0, 1).

[Jan Blechta]

Sure. Made a mistake. Only +-1 are nontrivial components of the gradient.

[Chris Richardson]

The 2017.1.0 release of ffc produces the same result (before merging Martin's changes, I think).

[Chris Richardson]

Sorry - my mistake - seems to be a new function, not in release...

[Prof Garth Wells]

It looks like the loop:

        for (std::size_t r = 0; r < num_derivatives; ++r)
        {
            // Declare derivative matrix (of polynomial basis).
            double dmats[3][3];
            // Initialize dmats.
            std::size_t comb = combinations[order - 1][r][0];
            for (std::size_t t = 0; t < 3; ++t)
                for (std::size_t u = 0; u < 3; ++u)
                    dmats[t][u] = dmats0[comb][t][u];
            // Looping derivative order to generate dmats.

needs to be brought inside the loop:

            case 0:
                // Compute reference derivatives for dof 0.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {

I.e.,

            case 0:
                // Compute reference derivatives for dof 0.
                for (std::size_t r = 0; r < num_derivatives; ++r)
                {

            // Declare derivative matrix (of polynomial basis).
            double dmats[3][3];
            // Initialize dmats.
            std::size_t comb = combinations[order - 1][r][0];
            for (std::size_t t = 0; t < 3; ++t)
                for (std::size_t u = 0; u < 3; ++u)
                    dmats[t][u] = dmats0[comb][t][u];
            // Looping derivative order to generate dmats.