Repository: kinnala/scikit-fem

Issue:
FacetBasis on a set of facets with ElementComposite not working
```python
from skfem import *

mesh = MeshTri().refined(3)

facet_basis = FacetBasis(mesh, ElementTriP1() * ElementTriP1(), facets='left')
facet_basis.interpolate(facet_basis.zeros())
```
results in
```
ValueError: operands could not be broadcast together with shapes (32,3) (8,3) (32,3) 
```

Discussion:

Patch:
diff --git a/skfem/assembly/basis/abstract_basis.py b/skfem/assembly/basis/abstract_basis.py
index 63290000..d352694e 100644
--- a/skfem/assembly/basis/abstract_basis.py
+++ b/skfem/assembly/basis/abstract_basis.py
@@ -295,13 +295,13 @@ class AbstractBasis:
 
         # loop over solution components
         for c in range(len(refs)):
-            ref = refs[c][1].basis[0][0]
-            ref = ref.astuple
-            fs = []
+            ref = refs[c][1].basis[0][0].astuple
 
-            def linear_combination(n, refn):
+            def linear_combination(n):
                 """Global discrete function at quadrature points."""
-                out = np.zeros_like(refn, dtype=w.dtype)
+                out = 0. * np.einsum('...,...j->...j',
+                                     w[self.element_dofs[0]],
+                                     self.basis[0][c].get(n))
                 for i in range(self.Nbfun):
                     values = w[self.element_dofs[i]]
                     out += np.einsum('...,...j->...j', values,
@@ -309,9 +309,10 @@ class AbstractBasis:
                 return out
 
             # interpolate DiscreteField
+            fs = []
             for n in range(len(ref)):
                 if ref[n] is not None:
-                    fs.append(linear_combination(n, ref[n]))
+                    fs.append(linear_combination(n))
                 else:
                     fs.append(None)
 


Test patch:
diff --git a/tests/test_basis.py b/tests/test_basis.py
index d77533ad..c39fd890 100644
--- a/tests/test_basis.py
+++ b/tests/test_basis.py
@@ -117,6 +117,13 @@ class TestCompositeFacetAssembly(TestCase):
         assert_allclose(A[0].todense()[0, ::2],
                         C[0].todense()[0])
 
+        y = fbasis1.zeros()
+        y[fbasis1.get_dofs('left')] = 1
+        assert_allclose(
+            fbasis1.project(fbasis1.interpolate(fbasis1.ones())),
+            y,
+        )
+
 
 class TestFacetExpansion(TestCase):