Test fiat #404

[ReubenHill]

No description provided.

[dham]

Hell yes. This is much better.

[ReubenHill]

I'll stick it in its own PR

[ReubenHill]

#2741

[ReubenHill]

I test the following

from firedrake import *
m = ExtrudedMesh(UnitSquareMesh(2, 3), 2)
V = FunctionSpace(m, "CG", 2)
f = Function(V)
f.at([0.5, 0.5, 0.4])

With fiat master, our kernel has

static inline void to_reference_coords_kernel(void *result_, double *x0, int *return_value, double *C)
{
    ...
    // Are we inside the reference element?
    *return_value = PetscRealPart(X[0]) + 1.0e-14 >= 0 && PetscRealPart(X[1]) + 1.0e-14 >= 0 && PetscRealPart(X[2]) + 1.0e-14 >= 0 && PetscRealPart(X[2]) - 1.0e-14 <= 1 && PetscRealPart(X[0]) + PetscRealPart(X[1]) - 1.0e-14 <= 1;
}

(16 operations, all adds and logical)

With branch [ReubenHill/point-locate](https://github.com/firedrakeproject/fiat/tree/ReubenHill/point-locate) (firedrakeproject/fiat#33) we now have

static inline void to_reference_coords_kernel(void *result_, double *x0, int *return_value, double *C)
{
    ...
    // Are we inside the reference element?
    *return_value = 0.5*fabs(PetscRealPart(X[2])) + 0.5*fabs(PetscRealPart(X[2]) - 1.0) - 0.5 <= 1.0e-14 && 0.5*fabs(PetscRealPart(X[0])) + 0.5*fabs(PetscRealPart(X[1])) + 0.5*fabs(PetscRealPart(X[0]) + PetscRealPart(X[1]) - 1.0) - 0.5 <= 1.0e-14;
}

(21 operations: 5 multiplications, 5 fabs, 8 add/sub, 3 logical)

[ReubenHill]

If we apply the following to [ReubenHill/point-locate](https://github.com/firedrakeproject/fiat/tree/ReubenHill/point-locate)/(firedrakeproject/fiat#33)

diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py
index fa4d47b..c4e880b 100644
--- a/FIAT/reference_element.py
+++ b/FIAT/reference_element.py
@@ -657,7 +657,7 @@ class UFCSimplex(Simplex):
         # sum(bary-abs(bary)) = -4.
         # - 0.5 * sum(bary-abs(bary)) = 2.0
         # which is the correct L1 distance from the cell to the point.
-        return - 0.5 * sum([b - abs(b) for b in bary])
+        return abs(0.5 * sum([b - abs(b) for b in bary]))
 
 
 class Point(Simplex):

(i.e. swap - for abs) we get

static inline void to_reference_coords_kernel(void *result_, double *x0, int *return_value, double *C)
{
    ...
    // Are we inside the reference element?
    *return_value = fabs(0.5*fabs(PetscRealPart(X[2])) + 0.5*fabs(PetscRealPart(X[2]) - 1.0) - 0.5) <= 1.0e-14 && fabs(0.5*fabs(PetscRealPart(X[0])) + 0.5*fabs(PetscRealPart(X[1])) + 0.5*fabs(PetscRealPart(X[0]) + PetscRealPart(X[1]) - 1.0) - 0.5) <= 1.0e-14;
}

(23 operations: 5 multiplications, 7 fabs, 8 add/sub, 3 logicals) which is slightly worse