-
Notifications
You must be signed in to change notification settings - Fork 159
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Test fiat #404 #2738
Test fiat #404 #2738
Conversation
@@ -65,8 +66,7 @@ def is_affine(ufl_element): | |||
def inside_check(fiat_cell, eps, X="X"): | |||
dim = fiat_cell.get_spatial_dimension() | |||
point = tuple(sympy.Symbol("PetscRealPart(%s[%d])" % (X, i)) for i in range(dim)) | |||
|
|||
return " && ".join("(%s)" % arg for arg in fiat_cell.contains_point(point, epsilon=eps).args) | |||
return ccode(fiat_cell.contains_point(point, epsilon=eps)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Hell yes. This is much better.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I'll stick it in its own PR
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
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 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) |
If we apply the following to 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 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 |
No description provided.