Attempt 1 at replicated Forward Mag Operator

ForwardMag.m

%% Geometry
nx = 9;
ny = nx;
nz = nx-1;

% Set up graded mesh
curve = 0;
weight = 0.61; % Percent finer mesh at finest scale to coarsest.
gmesh = meshfunc(curve, weight);

% domain limit
domain = 1;

% xyz on the cell nodes
[xn, hx] = gmesh(linspace(-domain, domain, nx));
[yn, hy] = gmesh(linspace(-domain, domain, ny));
[zn, hz] = gmesh(linspace(-domain, domain, nz));

% Helper vars
ncells = nx * ny * nz;
nfaces = (nx+1) * (ny) * (nz) + (nx) * (ny+1) * (nz) + (nx) * (ny) * (nz+1);
vol = kron(hz, kron(hy, hx));

% number of data
nobs = 10;

% Operators
G = getop(hx, hy, hz, 'gradc2f', 'dirichlet'); % Centres -> Faces
Av = getop(hx, hy, hz, 'avgfaces2centres', 'vector');
Ac = getop(hx, hy, hz, 'avgfaces2centres', 'scalar');

% Random independent variables
k = abs(randn(ncells,1));
h = abs(randn(nfaces,1));
R = abs(randn(nobs*ncells, 1));

% Base matrices
I = speye(ncells);
E = [I,I,I];
e = ones(nobs,1);

% Kron up operators to Nobs dimension
Gb = kron(speye(nobs), G);
Avb = kron(speye(nobs), Av);
Acb = kron(speye(nobs), Ac);
Vol = kron(speye(nobs), vol');

% Forward Op
Q = @(k) kron(e, E * (Av * (G*k) .* (Av*h)) ) ./ R;

F = @(k) Vol * (Acb * Gb * Q(k));

% Get some data
data = F(k);

Explanation

The operators are kron'd out in large block identities. So that Q(k) produces an nobs*ncell length vector. Q(k) is doing the inner integral part. Then it gets hit by the big Gb but that puts the vectors on the walls again. But since this is TMI or whatever, I Acb them, average them to the centre.... maybe I don't want to average, maybe just sum? And then I just multiply by Vol ...

The problem might be that Vol needs to be applied before averaging or summing the final Gradient.

Also I have no idea yet how to differentiate this sucker.