Fix restrict=True for actions

[leo-collins]

We need to replace the contracted argument

[connorjward]

Seems to make sense.

[pbrubeck]

What error were you getting? It should be possible still to assemble the original unrestricted F and then interpolate the result into the restricted space. This seems like a workaround

[leo-collins]

What error were you getting?

An error in UFL for incompatible spaces in action.

It should be possible still to assemble the original unrestricted F and then interpolate the result into the restricted space.

I don't understand. Where are we assembling F?

[pbrubeck]

I don't understand. Where are we assembling F?

assemble(interpolate(v_res, F(u->u_res; v)))

will break into assemble(F(u->u_res; v)) into a Cofunction(V.dual()), and then it will interpolate it to V_res.dual(). There is no need to change the internal arguments in F(u_res; v).

My guess is that the error in UFL might be doing the checks wrong.

[pbrubeck]

Also another potential source of issues is the fact that action with an interpolate are eagerly simplifying, maybe that is the problem.

[leo-collins]

We have F = action(action(M, I), u) - L (looking at the test I added). If we replace(F, {self.u: self.u_restrict}) then we end up with action(action(M, I), u_restricted). action(M, I) goes from V -> V^*, and u_restricted lives in V_restricted. Hence the incompatible spaces.

In compute_bc_lifting we do F = ufl_expr.action(J, u), where J=action(M, I). However this returns an Action(Action(M, I), u) since J is a BaseForm. If J is a regular Form, then ufl_expr.action(J,u) eagerly computes the action and returns another Form.

[pbrubeck]

We have F = action(action(M, I), u) - L (looking at the test I added). If we replace(F, {self.u: self.u_restrict}) then we end up with action(action(M, I), u_restricted). action(M, I) goes from V -> V^*, and u_restricted lives in V_restricted. Hence the incompatible spaces.

In compute_bc_lifting we do F = ufl_expr.action(J, u), where J=action(M, I). However this returns an Action(Action(M, I), u) since J is a BaseForm. If J is a regular Form, then ufl_expr.action(J,u) eagerly computes the action and returns another Form.

Sounds like one mathematically safe fix is to replace u -> interpolate(u_res, V)

[pbrubeck]

Suggested change

	self.F = interpolate(v_res, replace(F, {u_arg: u_res, self.u: self.u_restrict}))
	self.F = interpolate(v_res, replace(F, {self.u: interpolate(self.u_restrict, V)}))

[pbrubeck]

There are many ways to "fix" this. Of course my suggestion is far from optimal because we need to generate code and incur extra flops. But I think this is the most robust solution. We can make it efficient if we intercept V_res->V in interpolate and replace the codegen with dof-permutation (i.e. assign)

[pbrubeck]

Should we turn this into a reusable function? It might be good to generalize to any number of arguments