I constructed a MWE below for an example of the harmonic heat flow into the sphere with the weak formulation

(\partial_t d, \phi) - (d \times \Delta d, d\times \phi) =0

Hereby, I compute one iteration of three different variational formulations.

As matrix we simply consider the cross product in this example.

**Case 1** corresponds to the standard L^2 formulation.

(\partial_t d, \phi)_h - (d \times \Delta d, d\times \phi)_2 =0

with

(f,g)_h = \int_\Omega \mathcal{I} ( f\cdot g) dx

This corresponds to

```
inner(d - d0, td)*dxL - dt*inner(cross(d0,q), cross(d0,td)) *dx
```

**Case 3** corresponds to the case with mass-lumping applied.

(\partial_t d, \phi)_h - (d \times \Delta d, d\times \phi)_h =0

```
inner(d - d0, td)*dxL - dt*inner(cross(d0,q) , cross(d0,td)) *dxL
```

**Case 2** corresponds to the above described need of interpolating the Test Function with a matrix.

(\partial_t d, \phi)_h - (\mathcal{I} [d \times \Delta d], \mathcal{I} [d\times \phi ])_2 =0

I try to achieve this using the ufl Function `as_vector`

. Since this is part of the ufl standard repertoire I however expect the same behaviour as in case 1 due to the assembly.

MWE:

```
import numpy as np
from dolfinx.fem import Function, functionspace, form, assemble_scalar, ElementMetaData
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import CellType, create_box, CellType
from dolfinx.io import XDMFFile
from ufl import dx, grad, inner, div, dot, TrialFunctions, TrialFunction, TestFunction, TestFunctions, lhs, rhs, Measure, as_vector, cross, VectorElement, MixedElement, split
from mpi4py import MPI
def L2error(a):
error = form(inner(a, a) * dx)
error_a = assemble_scalar(error)
return error_a
def ic(x: np.ndarray)-> np.ndarray:
# x hase shape (dimension, points)
dim = 3
values = np.zeros((dim, x.shape[1])) # values is going to be the output
values[0]= np.sin( 2.0*np.pi*(np.cos(x[0])-np.sin(x[1]) ) )
values[1]= np.cos( 2.0*np.pi*(np.cos(x[0])-np.sin(x[1]) ) )
if dim == 3: values[2]=0.0
return values
def main():
dim = 3
dt = 0.5
dh = 10
domain = create_box(MPI.COMM_WORLD, [np.array([-1, -1,-1]), np.array([1.0, 1.0,1.0])], [dh,dh,dh], cell_type = CellType.tetrahedron)
el = VectorElement("Lagrange", domain.ufl_cell(), 1, dim = dim)
me = MixedElement([el,el])
FS = functionspace(domain, me)
d, q = TrialFunctions(FS)
td, tq = TestFunctions(FS)
u1 = Function(FS)
u2 = Function(FS)
u3 = Function(FS)
u0 = Function(FS)
d0, q0 = split(u0)
def g_times_f(g,f):
"""
cross product as matrix
https://en.wikipedia.org/wiki/Cross_product#Alternative_ways_to_compute
"""
f = as_vector([
-g[2]*f[1] + g[1]*f[2],
g[2]*f[0] - g[0]*f[2],
-g[1]*f[0]+g[0]*f[1]
])
return f
dxL = Measure("dx", metadata = {"quadrature_rule": "vertex"}) # mass lumping
# Example:
# Harmonic Heat Flow into sphere
energy_eq = inner( grad(d), grad(tq))*dx - inner(q,tq)*dxL
# Case 1: L2 case (standard)
eq1 = inner(d - d0, td)*dxL - dt*inner(cross(d0,q), cross(d0,td)) *dx + energy_eq
# Case 2: reformulation with ufl.as_vector
# Expected Behaviour: Same as L2 case 1
eq2 = inner(d - d0, td)*dxL - dt*inner(g_times_f(d0,q) , g_times_f(d0,td)) *dx + energy_eq
# Case 3: fully mass-lumped
eq3 = inner(d - d0, td)*dxL - dt*inner(cross(d0,q) , cross(d0,td)) *dxL + energy_eq
a1=lhs(eq1)
b1=rhs(eq1)
a2=lhs(eq2)
b2=rhs(eq2)
a3=lhs(eq3)
b3=rhs(eq3)
u0.sub(0).interpolate(ic)
solver_param = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type": "mumps"}
LP1 = LinearProblem(a1,b1 , bcs=[], u=u1, petsc_options=solver_param)
LP2 = LinearProblem(a2,b2 , bcs=[], u=u2, petsc_options=solver_param)
LP3 = LinearProblem(a3,b3 , bcs=[], u=u3, petsc_options=solver_param)
LP1.solve()
LP2.solve()
LP3.solve()
print("diff between 1 and 2", L2error(u1.sub(0).collapse()-u2.sub(0).collapse()))
print("diff between 1 and 3", L2error(u1.sub(0).collapse()-u3.sub(0).collapse()))
print("diff between 2 and 3", L2error(u3.sub(0).collapse()-u2.sub(0).collapse()))
if __name__ == '__main__':
main()
```

The output is:

```
diff between 1 and 2 0.0
diff between 1 and 3 3.8015563928781586
diff between 2 and 3 3.8015563928781586
```

This shows that using the method `as_vector`

is assembled in the same way as the standard ufl formulation. This was also my expected behaviour.

So I am still wondering if one can construct a Trial Function where the nodal coefficients can be scaled dependent on another function. So for a Test/Trial Function `u`

that can be represented as

u = \sum_z u_z \psi_z,

I instead want to use

\mathcal{I}(Au) = \sum_z (A_z u_z) \psi_z,

as Test/Trial Function in the variational formulation.

I hope that can explain my question and the motivation behind it better.