ValueError: Shape Mismatch when Comparing Interpolated Solution to Numerical

1

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)
2

Based off /qa/13269/scalar-mixed-function-space-problem/, the solution was to set deepcopy=True via u.split(True).

1 Like