I am solving a set of coupled Maxwell’s equations and got the system working. To calculate the error, I use analytical expressions for the two fields and am trying to compare those against my numerical results at degrees of freedom (d.o.f.).

I am trying to evaluate at the d.o.f. because I am using Nedelec elements of the 1st kind, and so think that I cannot the following function (which works perfectly) does not give the complete picture:

```
compute_vertex_values()
```

However, it turns out that I get the following error when I try to compute the error. I do not understand why, because I am using function_space from the solution to define the space I want to interpolate in.

```
ValueError: operands could not be broadcast together with shapes (1115,) (2230,)
```

It appears as if the interpolated data is half that of the numerical solution every time. Could someone help me see why? I have seen a few posts on here that discuss this issue, but none of them seem to work. For example, */qa/9444/interpolation-on-mixed-functionspace-shapes-not-matching/?show=9447#a9447* suggests using *FunctionAssigner()* instead of *split()*, but my implementation of the attempt did not change anything (perhaps I did it wrong). Here is another related post: */qa/14036/initial-condition-in-mixed-function-spaces/*.

Here is a minimal working example:

```
from dolfin import *
import numpy as np
[f, pie] = (1, np.pi)
[omega, k] = (Constant(2*pie*f), 2)
[eps, mu] = (Constant(1.0), Constant(1.0))
[Nx, Ny, Nz] = (5, 5, 5)
mesh = UnitCubeMesh(Nx, Ny, Nz)
P1 = FiniteElement("N1E", tetrahedron, 1) # FiniteElement("N1curl", mesh.ufl_cell(), 1) nedelec edge elements
TH = P1*P1
V = FunctionSpace(mesh, TH)
# Define test functions
test_E, test_H = TestFunctions(V)
# Define functions for velocity and concentrations
u = Function(V)
# Split system functions to access components
trial_E, trial_H = split(u)
F = omega * inner(eps*trial_E, test_H) * dx + omega * inner(mu * trial_E, test_E) * dx
solve(F == 0, u)#, bcs) # solve(a == L, u, bc). NO bc rn.
_u_1, _u_2 = u.split()
# Analytical solution:
EE = Expression(("cos(2*piee*kk*x[2])", "cos(2*piee*kk*x[2])", "0"), piee = pie, kk = k, degree=2)
HH = Expression(("-cos(2*piee*kk*x[2])", "cos(2*piee*kk*x[2])", "0"), piee = pie, kk = k, degree=2)
# Interpolate analytical result:
V1 = _u_1.function_space().collapse()
V2 = _u_2.function_space().collapse()
exact_E = interpolate(EE, V1)
exact_H = interpolate(HH, V2)
nodal_values_EE = exact_E.vector().get_local()
nodal_values_HH = exact_H.vector().get_local()
# Compute error
error_max_EE = (nodal_values_EE - _u_1.vector().get_local()).max()
error_max_HH = (nodal_values_HH - _u_2.vector().get_local()).max()
print('error_max_EE =', error_max_EE)
print('error_max_HH =', error_max_HH)
```