Certified touching-panel and disjoint-panel quadrature for the Laplace P0 matrix

Parametrize $T$ and $S$ by affine maps $t(u)$ and $s(v)$ from ${\mathrm{\Delta }}_2=\{(u_1,u_2):u_i\ge 0,u_1+u_2\le 1\}$. Their constant surface Jacobians are $J_T=\mathrm{\mid }(t_1-t_0)\times (t_2-t_0)\mathrm{\mid }$ and $J_S=\mathrm{\mid }(s_1-s_0)\times (s_2-s_0)\mathrm{\mid }$. With $L(u,v)=t(u)-s(v)$ on $P={\mathrm{\Delta }}_2\times {\mathrm{\Delta }}_2$, the nonnegative change of variables in [F8] gives ${\mathcal{I}}_{TS}=(J_T{J}_S\mathrm{/}(4\pi )){\int }_P\mathrm{\mid }L(z){\mathrm{\mid }}^{-1}dz$, initially allowing an infinite value. The zero set $Z=P\cap L^{-1}(0)$ is the parameterized diagonal of $T\cap S$ and has dimension $2$, $1$, or $0$ for coincident, edge-sharing, or vertex-sharing panels.

Cut $P$ by each non-identically-zero coordinate hyperplane $L_i=0$. Enumerate each cell’s finitely many vertices by active boundary-hyperplane sets, discarding infeasible intersections. The at most eight resulting full-dimensional sign cells cover $P$ up to their boundaries; those boundaries are contained in finitely many hyperplanes and have four-dimensional measure zero by Fubini in [F8]. For any $z_0\in Z$ and interior point $z$ of a sign cell, $L((1-t)z_0+tz)=tL(z)$ for $0<t\le 1$, so every $z_0$ belongs to that cell's closure. The simultaneous zero set of its signed nonnegative affine coordinates is therefore a nonempty face.

Triangulate each cell by pulling from a global vertex ordering that lists all vertices in $Z$ first and uses the same induced order on shared faces. The pulling induction is complete: for any point other than the first vertex, extend the segment from that vertex through the point to its last point in the cell; the exit point lies on a facet not containing the first vertex, so cones over those facets cover the cell, and overlaps occur only on their boundaries. Induction triangulates those facets and the common induced ordering makes adjacent cell triangulations agree. Every full-dimensional cell has a vertex in its nonempty face $Z$, hence its first vertex is in $Z$; each resulting four-simplex therefore has a zero vertex. In a simplex with signed affine coordinates nonnegative, their common zero set is exactly the convex hull of the zero vertices, so the zero face has dimension $j\in \{0,1,2\}$, covering all three contact types including coincident panels.

In a resulting four-simplex $\mathrm{\Sigma }=[v_0,\dots ,v_4]$, order the vertices so exactly $v_0,\dots ,v_j$ lie in $Z$, put $d=3-j$ and $D=\mathrm{\mid }\det (v_1-v_0,\dots ,v_4-v_0)\mathrm{\mid }$, and define $B(b)={\sum }_{m=j+1}^4{b}_{m}L(v_m)$ on the opposite face ${\mathrm{\Delta }}_d$. That face contains no zero of $L$, so compactness gives $c={\min }_{{\mathrm{\Delta }}_d}\mathrm{\mid }B\mathrm{\mid }>0$. The join coordinates $z=(1-r){\sum }_{i=0}^j{a}_i{v}_i+r{\sum }_{m=j+1}^4{b}_m{v}_m$, with $a\in {\mathrm{\Delta }}_j$, $b\in {\mathrm{\Delta }}_d$, and $0\le r\le 1$, are one-to-one in the interiors and have absolute Jacobian $D{r}^{3-j}(1-r{)}^j$: the $j$ independent $a$ directions scale by $1-r$, the $d$ independent $b$ directions by $r$, and the radial column supplies $D$. Also $L(z)=rB(b)$.

The kernel times this Jacobian has radial factor $r^{2-j}(1-r{)}^j$. For $j=0,1,2$, ${\int }_0^1{r}^{2-j}(1-r{)}^{j}dr=(2-j)!j!\mathrm{/}3!$ and $\mathrm{\mid }{\mathrm{\Delta }}_j\mathrm{\mid }=1\mathrm{/}j!$, so ${\int }_{\mathrm{\Sigma }}\mathrm{\mid }L(z){\mathrm{\mid }}^{-1}dz=D(2-j)!\mathrm{/}6{\int }_{{\mathrm{\Delta }}_d}\mathrm{\mid }B(b){\mathrm{\mid }}^{-1}db$. This is finite because $\mathrm{\mid }B\mathrm{\mid }\ge c$. For $d=1$, split the unit square into its two triangles and swap coordinates on one to get ${\int }_0^{1}g(t)dt={\int }_R(g(x)+g(y))dxdy$ with $g=\mathrm{\mid }B{\mathrm{\mid }}^{-1}$. For $d=2$, $(u,v)=(1-x,y)$ maps $R=\{0\le y\le x\le 1\}$ onto ${\mathrm{\Delta }}_2$ with determinant one. For $d=3$, set $u=1-x$, $v=y$, $w=(x-y)t$ on $R\times [0,1]$; the remaining barycentric coordinate is $(x-y)(1-t)$ and the Jacobian is $x-y$. Thus each case is a finite sum of reference-triangle integrals, including all zero-face dimensions, and each transformed amplitude is smooth near $R$ because its denominator is bounded below by $c$. The coefficients $c_m$ include the positive factor $J_T{J}_S\mathrm{/}(4\pi )$ from step 1.1 and the positive simplex determinants and radial factors above.

The positive bound $c$ is effectively certifiable from the finitely many vectors $L(v_m)$. Enumerate every nonempty face of their convex hull and solve the finite affine normal equations for the point nearest zero on its affine span, retaining feasible barycentric coordinates; a global minimizer lies in the relative interior of a smallest face and has zero tangential derivative there, so this enumeration finds the minimum. Row-reduce degenerate systems and include their lower-dimensional faces. Positivity follows from the zero-face characterization in step 4.1. For algebraic panel coordinates, exact algebraic comparisons and interval refinement yield a rational $0<\underset{‾}{c}<c$. For each polynomial parameterization $B(\xi )$ whose admissible domain lies in a unit cube and whose image lies in the image under $L$ of the opposite face, set $h=\mathrm{\mid }B{\mathrm{\mid }}^2\ge c^2$ on that admissible domain, $P_{0,0}=1$, out-of-range $P_{\beta ,j}=0$, and $P_{\beta +e_i,j}={\mathrm{\partial }}_i{P}_{\beta ,j}-(j-\frac{1}{2})P_{\beta ,j-1}{\mathrm{\partial }}_{i}h$. The product and chain rules inductively give ${\mathrm{\partial }}^{\beta }{h}^{-1\mathrm{/}2}={\sum }_{j=0}^{\mathrm{\mid }\beta \mathrm{\mid }}{P}_{\beta ,j}{h}^{-1\mathrm{/}2-j}$; hence $M_{\beta }={\sum }_j\mathrm{coefabs}(P_{\beta ,j}){\underset{‾}{c}}^{-1-2j}$ bounds each required derivative there, since the absolute coefficient sum bounds a polynomial on the containing unit cube. Compose the triangle coordinates with the Duffy map before using a full square as the admissible domain; no positive denominator bound is asserted on an extension outside that image. Apply the finite Leibniz rule to polynomial prefactors such as $x-y$ and add the bounds for finite amplitude sums before using these estimates for $q_m$. Differentiation under the finite $t$ integral in the $d=3$ map is justified by these uniform bounds and the fundamental theorem of calculus.

For each smooth reference amplitude $q_m$ in step 5.1, the Duffy map $D(r,s)=(r,rs)$ from [F5] gives ${\int }_R{q}_m={\int }_0^1{\int }_0^{1}r{G}_m(r,s)dsdr$, $G_m(r,s)=q_m(r,rs)$. With $n\ge 1$, the positive product rule in [F6] has error at most $K_1\mathrm{\parallel }{\mathrm{\partial }}_r^{2n}{G}_m{\mathrm{\parallel }}_{\mathrm{\infty }}+\frac{1}{2}{K}_0\mathrm{\parallel }{\mathrm{\partial }}_s^{2n}{G}_m{\mathrm{\parallel }}_{\mathrm{\infty }}$, since ${\int }_0^{1}rdr=1\mathrm{/}2$. The derivative bounds are computed in step 6.1; add the separately certified node-value, fiber, and outward-arithmetic radii with the coefficients $c_m$.

This bound gives effective termination. Fix $n=1$ and subdivide $R$ into similar triangles of maximum diameter $h$; after affine pullback the second pure derivatives of $q_m(A(r,rs))$ are bounded by $C{h}^2\mathrm{\parallel }{D}^2{q}_m{\mathrm{\parallel }}_{\mathrm{\infty }}$, so the positive-rule errors sum to $O(h^2)$ after multiplying by the affine Jacobians, whose sum is one. Indeed, $D(r,s)=(r,rs)$ has ${\mathrm{\partial }}_r^{2}D={\mathrm{\partial }}_s^{2}D=0$, so the pure second-derivative chain rules contain only the Hessian of $q_m$ and two factors from the affine map, each of size $O(h)$. Its nonzero mixed derivative is not part of the Gaussian remainder used here. In the $d=3$ fiber, composite midpoint quadrature on $M$ equal subintervals has error at most $B_{tt}\mathrm{/}(24M^2)$ by the one-variable Gaussian remainder [F7]. Choose $h,M$ so these bounds consume less than half the requested tolerance, then refine outward interval enclosures for the remaining half. Exact algebraic input makes incidence, determinant, and strict-positive denominator decisions finite; every refinement terminates because the certified lower bounds are positive. With other coordinate formats, require supplied certified incidence, nondegeneracy, separation, and arithmetic bounds; absent or failed certificates return uncertified.

If $T,S$ are disjoint, $d_0=\mathrm{dist}(T,S)>0$. The parameter kernel $g(u,v)=\mathrm{\mid }L(u,v){\mathrm{\mid }}^{-1}$ is smooth on ${\mathrm{\Delta }}_2\times {\mathrm{\Delta }}_2$ with denominator at least $d_0$. Compose $L$ with the Duffy parameterizations of both triangles; it is polynomial on a product of unit squares and still has norm at least $d_0$, so the recurrence of step 6.1 with a rational lower bound below $d_0$ supplies explicit uniform derivatives. Integrating over $v$ defines $q(u)={\int }_{{\mathrm{\Delta }}_2}g(u,v)dv$; differentiation under this finite integral is justified by those uniform bounds, and the simplex mass is $1\mathrm{/}2$. The affine map $(u_1,u_2)=(1-x,y)$ sends the outer triangle to $R$. Evaluate each inner $v$ integral by the positive composite $\alpha =0$ Duffy rule; its error bound is uniform in $u$. Apply the outer rule of step 7.1 to $q$. Subdividing both reference triangles and using their computed second-derivative bounds makes the inner errors $O(h^2)$ uniformly in $u$ and the outer error $O(h^2)$; positivity and bounded total weight let these errors add. Thus the disjoint smooth fallback terminates at every positive tolerance, with its two numerical radii charged separately.

To identify the double integral with the weak matrix entry, use barycentric coordinates $(u_1,u_2,1-u_1-u_2)$ and cutoffs ${\chi }_n={\prod }_{i=1}^3\min (1,n{u}_i)$ on each panel. They are bounded by one and converge almost everywhere and in $L^2$ to the panel indicator. By [F1], each cutoff is an $L^2$ density and the fixed $H^{1\mathrm{/}2}$ chart norm controls $L^2$, so these cutoffs converge in $H^{-1\mathrm{/}2}$ by duality; boundedness of $V$ in [F3] gives $⟨{\chi }_n{1}_T,V({\chi }_n{1}_S)⟩\to ⟨1_T,V{1}_S⟩$. For the trace identification, a finite union of planar faces satisfies $\mathrm{\mid }\mathrm{\Gamma }\cap B(x,r)\mathrm{\mid }\le C{r}^2$; dyadic annuli imply ${\int }_{\mathrm{\Gamma }\cap B(x,\delta )}\mathrm{\mid }x-y{\mathrm{\mid }}^{-1}d{S}_y\le C\delta$. Splitting at distance $2\delta$ and using $\mathrm{\mid }{\mathrm{\nabla }}_x{G}_3\mathrm{\mid }=(4\pi \mathrm{\mid }x-y{\mathrm{\mid }}^2{)}^{-1}$ gives continuity modulus $C\delta (1+\mathrm{\mid }\log \delta \mathrm{\mid })$ for each bounded cutoff potential throughout a bounded neighborhood of $\mathrm{\Gamma }$: the area and near-integral bounds are uniform for ambient target points, and the same near/far split applies there. By [F3] the off-surface kernel integral represents its weak potential; by [F4] its continuous restriction is its Sobolev trace, and Fubini [F8] identifies the pairing with the cutoff double integral. Absolute integrability from the touching-panel calculation in step 5.1 and the disjoint case in step 8.2 dominates these integrands, so dominated convergence gives $⟨1_T,V{1}_S⟩={\mathcal{I}}_{TS}$ for coincident, edge-sharing, vertex-sharing, and disjoint panels.