NotImplementedError when projecting during asm

[gatling-nrl]

I'm trying to work with something from P0 and something from P1 in the same form. (Specifically, in the form I want to do grad(w['thing_p1']) * w['thing_p0'] as part of a possible solution to #821.)

But I've run into a NotImplementedError and I'm not too sure what isn't implemented.

import skfem

mesh = skfem.MeshTri()

basis_p1 = skfem.Basis(mesh, skfem.ElementTriP1())
basis_p0 = basis_p1.with_element(skfem.ElementTriP0())

f1 = basis_p1.zeros()
f0 = basis_p0.zeros()

@skfem.Functional
def issue_demo(w):
    return w.x[0] # arbitrary, just something of the correct shape

skfem.asm(issue_demo, basis_p1, f1=f1, f0=basis_p1.project(f0))

raises

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
~/ipykernel_21860/2766063479.py in <module>
     13     return w.x[0] # arbitrary, just something of the correct shape
     14 
---> 15 skfem.asm(issue_demo, basis_p1, f1=f1, f0=basis_p1.project(f0))

\scikit-fem\skfem\assembly\basis\cell_basis.py in project(self, interp, elements)
    229         from skfem.utils import solve, condense
    230 
--> 231         M, f = self._projection(interp)
    232 
    233         if elements is not None:

\scikit-fem\skfem\assembly\basis\abstract_basis.py in _projection(self, interp)
    432         return (
    433             BilinearForm(lambda u, v, _: inner(u, v)).assemble(self),
--> 434             LinearForm(lambda v, w: inner(interp, v)).assemble(self),
    435         )
    436 

\scikit-fem\skfem\assembly\form\form.py in assemble(self, *args, **kwargs)
     77         assert self.form is not None
     78         logger.info("Assembling '{}'.".format(self.form.__name__))
---> 79         out = (COOData(*self._assemble(*args, **kwargs))  # type: ignore
     80                .todefault())
     81         logger.info("Assembling finished.")

\scikit-fem\skfem\assembly\form\linear_form.py in _assemble(self, ubasis, vbasis, **kwargs)
     42             ixs = slice(nt * i, nt * (i + 1))
     43             rows[ixs] = vbasis.element_dofs[i]
---> 44             data[ixs] = self._kernel(vbasis.basis[i], w, dx)
     45 
     46         return np.array([rows]), data, (vbasis.N,), (vbasis.Nbfun,)

\scikit-fem\skfem\assembly\form\linear_form.py in _kernel(self, v, w, dx)
     47 
     48     def _kernel(self, v, w, dx):
---> 49         return np.sum(self.form(*v, w) * dx, axis=1)

\scikit-fem\skfem\assembly\basis\abstract_basis.py in <lambda>(v, w)
    432         return (
    433             BilinearForm(lambda u, v, _: inner(u, v)).assemble(self),
--> 434             LinearForm(lambda v, w: inner(interp, v)).assemble(self),
    435         )
    436 

\scikit-fem\skfem\helpers.py in inner(u, v)
    109     elif len(U.shape) == 4:
    110         return ddot(U, V)
--> 111     raise NotImplementedError
    112 
    113 

NotImplementedError: 

[kinnala]

The error is very confusing, I'm sorry about that. This is the fix:

skfem.asm(issue_demo, basis_p1, f1=f1, f0=basis_p1.project(basis_p0.interpolate(f0)))

Moreover, we obviously need some checks in project so that this does not happen.

[kinnala]

The thing is that project assumes it needs to call basis_p1.interpolate when calling basis_p1.project with a numpy array even though f0 is defined using basis_p0.

[gatling-nrl]

The need for interpolate is confusing because I would have expected basis_p1.with_element to have resolved that. Is the underlying issue the fact that you can't tell what "kind" of numpy array you're working with? Is it worth considering subclassing ndarray (a non-trivial task to do correctly) so that it can carry some additional attributes, specifically its basis?

[kinnala]

I've been considering about the possibility of subclassing ndarray but have so far have decided against it. Working with the native objects of numpy and scipy is somehow fundamental to this library. All other libraries have some kind of GridFunction and BilinearOperator instead of simple ndarray and scipy.sparse. While it would enable many automagic things I'd prefer keeping it this way.

[gatling-nrl]

I think this issue and #847 are looking like pretty much the same thing now. Not sure if you can merge them... or maybe just close one of them.

I think it boils down to DiscreteField needing a little more to duck type ndarray and/or the question of whether to subclass ndarray for discrete fields and projections. (Do you think "discrete fields" and "projections" are the right words for the two kinds of data structures we're talking about here?)