Issue using a surfacefun as a pde coefficient

[scottweady]

There is an issue when using a surfacefun as a coefficient in a surfaceop. Here is a mwe:

p = 12;
nref = 2; %mesh should have more than one element
dom = surfacemesh.square(p + 1, nref);

% Works fine
c = 1;
pdo = struct('lap', 1, 'c', c);
L = surfaceop(dom, pdo);

% Throws an error
c = surfacefun(@(x, y, z) 0 * x + 1, dom);
pdo = struct('lap', 1, 'c', c);
L = surfaceop(dom, pdo);

The full error message is

Error using reshape
Number of elements must not change. Use [] as one of the size inputs to automatically calculate the appropriate size for that dimension.

Error in surfaceop.leaf.initialize_DtN (line 92)
        op.(name) = reshape([op.(name).vals{:}], [n^2 numPatches]);

Error in surfaceop/initialize (line 20)
        S.patches = surfaceop.leaf.initialize_DtN(S.op, S.domain, varargin{:});

Error in surfaceop (line 62)
            obj.initialize(rhs);

Error in mwe (line 13)
L = surfaceop(dom, pdo);

The problem comes from .vals{:} and seems to be a Matlab issue. For example, in the same block as the error (line 92 in surfaceop.leaf.initialize_DtN.m) these three lines give different results

  size([op.(name).vals{:}]) % returns [13 13]
  
  tmp = op.(name);
  size([tmp.vals{:}]) % returns [13 13]
  
  tmp = op.(name).vals;
  size([tmp{:}]) % returns [13 208]

A simple fix is to replace .vals{:} with .vec() in both line 92 in initialize_DtN.m and line 96 in initialize_ItT.m:

op.(name) = reshape(op.(name).vec(), [n^2 numPatches]);