Dear @Nathan_Langlet,
Sorry for the confusion in the last message. There is a lot of stuff here, let’s try to unpack…
I mentioned JAX because a Jacobian was mentioned; that’s why you need two operations that go beyond UFL:
- Inverse of a high-order tensor
- Since 1. goes beyond UFL, you need to compute the derivative of \tau, which is inverse.
JAX is able to do both, that’s why.
Supposing that \tau = M^{-1} and you consider then a term like v^TM^{-1}w.
During the FEniCS conference, I mentioned that instead of computing M^{-1} each iteration for every Gauss point, maybe it makes sense to solve a small linear problem My = w, where y = M^{-1}w and then you compute v^T y instead. In this case, y is an external operator indeed, but w is not an operand. As for the derivative, if Q is your main field variable with respect to which you take a directional derivative and y(Q) = M^{-1}(Q)w(Q), you need to understand what y^\prime(Q) looks like and what you need to compute via JAX/NumPy/DOLFINx… In particular, you need to compute w(Q), which is your StrongRes. Since it’s perfectly expressible with UFL, express it with UFL! Then use fem.Expression to compute its values wherever you want. Then you can compute the derivative of w(Q) with UFL and then compute its values via fem.Expression. BTW, this concerns the computation of M as well. Then you wrap all this with a Python function, something like:
def y_impl(Q: np.ndarray):
w_values = w_expr.eval(domain, cells)
M_values = M_expr.eval(domain, cells)
y_values = solve_linear(M_values, w_values) #any code inside
return y_values.reshape(-1)
Then for the derivative, it’s gonna be more complicated, but possible I think!
All this looks pretty complex, but the reason here is really to reduce the size of the external operator and the allocated memory.
What I’d suggest now is to implement the direct matrix inverse without the extra complexity I’ve just introduced, i.e. now \tau is an external operator. This should be much more straightforward to implement, and thus you will get an idea of how external operators work. You can come back to optimiziation like suggested above, later.
But before doing anything, I’d really recommend explicitly deriving the variational problem you are solving, and the following may be useful.
I recently dropped an update for the documentation: Some Notation for External Operators — External operators within FEniCSx/DOLFINx. This introduces a notation around directional derivatives and external operators. Since the use case of external operators is when the variational formulations are really complex, it may be challenging to understand how the external operators play their role. I see it as a context for a modern agentic LLM like Claude or Gemini. You give it two files: your formulation and this notation dolfinx-external-operator/doc/notes/notation.md at main · a-latyshev/dolfinx-external-operator · GitHub. Then it may guide you through your problem using our package. Never tested though, hope it will be helpful 
