Tensor composite Simpson force evaluator on a rectangle

Input real endpoints $a,b,c,d$, integer panel counts $n_x,n_y$, and an exact evaluation oracle that, when queried at a point, returns either a value in ${\mathbb{R}}^2$ or evaluation-failure. First, if $a\ge b$ or $c\ge d$, return invalid-domain without querying the oracle. Next, if either panel count is odd or is less than two, return invalid-panel-count without querying the oracle. For the remaining inputs set $R=[a,b]\times [c,d]$, $h_x=(b-a)\mathrm{/}{n}_x$, $h_y=(d-c)\mathrm{/}{n}_y$, $x_i=a+i{h}_x$ for $0\le i\le n_x$, and $y_j=c+j{h}_y$ for $0\le j\le n_y$. Query the oracle at every grid point in lexicographic order with $i$ first and $j$ second, storing the returned values. If a query fails, stop at the first such point and return evaluation-failure with that point. Otherwise, for each $j$ form the continuous piecewise-linear interpolant through the stored values $F(x_i,y_j)$ and apply Composite Simpson rule on an even uniform grid componentwise in $x$ to this interpolant, obtaining $G_j=Q_{x}F(\cdot ,y_j)$. Form the continuous piecewise-linear interpolant through the stored $G_j$ and apply the same valid Simpson procedure componentwise in $y$ to obtain $Q_{y}G=Q_y{Q}_{x}F$. Return this vector $Q_y{Q}_{x}F$. Each interpolant is queried only at its stored grid nodes, so no new oracle call is needed. Both calls to the supplier use strictly ordered endpoints and even panel counts at least two. This is an exact-real procedure; termination and force-accuracy bounds are established in Tensor composite Simpson error for a two-dimensional force integral ↗, and floating-point roundoff is outside its scope.