Product of two functions then projection on the functionspace leads to strange behavior

1

Hi all, I am wondering how does product of two function work in Fenics?

Followed is a snapshot of the code

P = FunctionSpace(mesh, "CG", 1) 
a = Function(P)
d2v_low = dof_to_vertex_map(P)
a.vector()[:] =init_gt_1[d2v_low]
b =  Function(P)
d2v_low = dof_to_vertex_map(P)
b.vector()[:] =init_gt_2[d2v_low]

evaluate_coords_1 = (a*b)(coords)
holder = project(a*b, P)
evaluate_coords_2 = holder(coords)

I thought evaluate_coords_1==evaluate_coords_2 but it is not the case. Any thoughts on why this might be happening?

2

a and b are CG-1 functions, i.e. piecewise linear. Their product is therefore CG-2, piecewise quadratic. Projecting onto P, a CG-1 space, finds the best piecewise linear fit to a piecewise quadratic function. This representation is clearly not exact. Consider the following MWE:

from fenics import *

mesh = UnitIntervalMesh(4)
P = FunctionSpace(mesh, "CG", 1) 
a = Function(P)
b = Function(P)

expr = Expression('x[0]', degree=1)
a.interpolate(expr)
b.interpolate(expr)

coords=[0.5]
evaluate_coords_1 = (a*b)(coords)

holder = project(a*b, P)
evaluate_coords_2a = holder(coords)

P2 = FunctionSpace(mesh, "CG", 2) 
holder = project(a*b, P2)
evaluate_coords_2b = holder(coords)

print(evaluate_coords_1)
print(evaluate_coords_2a)
print(evaluate_coords_2b)

0.25
0.23958333333333337
0.25000000000000006

You can see that projecting onto a CG-2 space gives much better accuracy.

2 Likes