Proof. The argument is the same for both formulations, so we focus on the former. First note that \(u^D_h|_E\) satisfies the required boundary condition on \(E\) (in a trivial manner). Consider the subset of nodal basis functions supported in \(\overline{E}\), that is, \(N_E = \{i\: :\: \mathrm{supp}(\varphi_i)\subseteq \overline{E}\}\). Since 14 holds for all \(v=\varphi_i\) with \(i\in N_E\), all the integrals in 14 are, in effect, restricted to \(E\). Also note that \(V_0^h(\overline{E})\) is generated by the set \(\{\varphi_i|_E\: :\: i\in N_E\}\). Hence, \(u^D_h|_E\) satisfies 14 with the set \(E\) replacing \(D\). ◻

Proof. Let \((f,u_h^b) \in (V_0^h(\overline{D}))^*\times V^h(\partial D)\) with \(f\) nonnegative, and \(u_h={\mathcal{S}}^{D}_h(f,u_h^b)\). Denote by \(m=\min_{i=1}^N u_h(P_i)\); this is the minimum among the values at the vertices, but for piecewise linear functions it coincides with the global minimum of \(u_h\) on \(\overline{D}\). Consider the set of vertices \[{\mathcal{C}}_m=\{Q\;\textrm{vertex\;in\;}\;{\mathcal{T}}_h\: : \: u_h(Q)=m\}.\] The arguments vary only slightly between (A), (B), and .

For (A), if \({\mathcal{C}}_m \subseteq \partial D\), then \(u_h\) attains its minimum only on \(\partial D\), and 17 holds in a strict sense. Otherwise, there exists \(Q\in{\mathcal{C}}_m \cap Int(D)\). It remains to show that \(u_h\) is constant, then the conclusion follows. Let \(E_{Q}\) be as in the hypothesis. Since the minimum of \(u_h\) on \(\overline{E}_{Q}\) is achieved at \(Q\) (because it is a global minimum of \(u_h\)), and by 16 \[\begin{align} u_h|_{E_Q} = {\mathcal{S}}^{E_Q}_h(f_{E_Q},u_h|_{\partial E_Q}), \end{align}\] it follows from property (A2) that \(u_h\) is constant on \(\overline{E}_Q\). Now let \(\widehat{Q}\in Int(D)\) be an arbitrary interior vertex, and let \(Q=Q_0, Q_1, \dots, Q_r=\widehat{Q}\) be a sequence of adjacent interior vertices connecting \(Q\) to \(\widehat{Q}\), that is, \(Q_{i-1}\) is adjacent to \(Q_{i}\) for \(i=1, \dots, r\). Since \(Q_i\in Int(E_{Q_i})\) and \(Q_{i-1}\) is adjacent to \(Q_i\), it follows that \(\{Q_{i-1},Q_i\} \subset E_{Q_{i-1}} \cap E_{Q_{i}}\). Assume \(u_h(Q_k)=m\) for some \(0\le k \le r-1\). Since \(m\) is the global minimum, the same argument used for \(k=0\) shows that \(u_h\) is constant on \(\overline{E}_{Q_{k}}\). Hence \(u_h(Q_{k+1})=m\). By induction over \(k\) we get \(u_h(Q_r)=u_h(Q_0)=m\). So \(u_h(\widehat{Q})=m\) for all interior vertices \(\widehat{Q}\). If \(R\) is a vertex lying on \(\partial D\), consider a vertex \(\widehat{Q}\in Int(D)\) that it is adjacent to \(R\). Since \(R\in E_{\widehat{Q}}\), it follows that \(u_h(R)=m\) as well, showing that \(u_h\) is constant.

For (B), if \({\mathcal{C}}_m \subseteq \partial D\), then \(u_h\) attains its minimum only on \(\partial D\). Hence, if \(P\) is an interior vertex, then \[\begin{align} u_h(P) > m = \min u_h^b = \min \left((u_h^b)^+-(u_h^b)^- \right) \ge \min -(u_h^b)^- = -\max (u_h^b)^-, \end{align}\] showing a strict inequality holds in 18 . Independently of the condition \({\mathcal{C}}_m \subseteq \partial D\), if \(m>0\), then \[\begin{align} u_h(P) \ge m >0 \ge \min -(u_h^b)^- = -\max (u_h^b)^-, \end{align}\] which shows, again, that 18 holds strictly. This leaves us with the case when \({\mathcal{C}}_m \cap Int(D) \ne \varnothing\) and \(m\le 0\). Let \(Q\in{\mathcal{C}}_m \cap Int(D)\). As in case (A), we will show that \(u_h\) is constant on the set \(\overline{E}_{Q}\) from the hypothesis, and the rest of the argument follows the same path as in the case (A). The focus is on \(u_h|_{E_Q}\), to which we apply the conditions (B1) and (B2). By (B1) we have \[\begin{align} \label{eq:m1ineq} 0\ge m = u_h(Q) \ge -\max (u_h^b)^- = m_1 = -(u_h^b)^-(Q_1), \end{align}\tag{20}\] for some \(Q_1\in \partial E_Q\). If \(m_1=0\), then by 20 we have \(m=m_1=0\); hence, Condition (B2) implies \(u_h|_{E_Q}\) is constant. If \(m_1<0\), then \((u_h^b)^-(Q_1) = -m_1 > 0\), which implies (by the definition of \((u_h^b)^-\)) that \((u_h^b)^-(Q_1) = -u_h^b(Q_1)\), showing that \(u_h^b(Q_1) = m_1\). Cf. 20 \[m = u_h(Q) \ge u_h(Q_1) = u_h^b(Q_1)= m_1.\] Since \(m\) is the global minimum of \(u_h\), it follows that \(m=m_1\). Now (B2) implies that \(u_h\) is constant on \(\overline{E}_Q\).
The proof for (R) is identical to that of (A), except we assume from the beginning that \(m\ge 0\). ◻

Proof. This follows easily from Theorem 4, since for each interior vertex \(Q\) there exists \(i\in I\) so that \(Q\in Int(E_i)\), and \({\mathcal{S}}^{E_i}_h\) satisfies sDMP-X. Conditions (A1)–(A2) (or (B1)–(B2), for \(X=B\)) are clearly verified, since they are weaker than those of sDMP-A and sDMP-R (or SDMP-B). ◻

Proof. We restrict our attention to the case when \(f\) in 14 is nonnegative. The matrix form of 14 with \(c\equiv 0\) is \[\label{eq:matforFEczero} {\mathbf{A}}_{\alpha \alpha} {\mathbf{u}}_{\alpha} + {\mathbf{A}}_{\alpha \beta} {\mathbf{u}}_{\beta} = {\mathbf{F}}_{\alpha}.\tag{24}\] Note that \({\mathbf{A}}{\mathbf{1}} = {\mathbf{0}}\), which implies \[\label{eq:stiffnessone} {\mathbf{A}}_{\alpha \alpha} {\mathbf{1}}_{\alpha} + {\mathbf{A}}_{\alpha \beta} {\mathbf{1}}_{\beta} = {\mathbf{0}}_{\alpha}.\tag{25}\] Hence, \[\begin{align} \label{eq:1alpha} {\mathbf{1}}_{\alpha} = -{\mathbf{A}}_{\alpha \alpha}^{-1}{\mathbf{A}}_{\alpha \beta} {\mathbf{1}}_{\beta}. \end{align}\tag{26}\] Assume 23 holds. Let \(u_h^b \in V^h(\partial D)\) be arbitrary, with vector representation \({\mathbf{u}}_{\beta}\in \mathbb{R}^{|\beta|}\), and \(u_h = {\mathcal{S}}^D_h(f,u_h^b)\), with vector representation (in the nodal basis) \({\mathbf{u}}\). Let \(\underline{u} = \min u_h^b = \min_i \{{\mathbf{u}}_i\;:\;i\in \beta\}\). Cf. 24 and 26 we have \[\begin{align} \nonumber {\mathbf{u}}_{\alpha}-\underline{u}{\mathbf{1}}_{\alpha}& = & {\mathbf{A}}_{\alpha \alpha}^{-1}\left( {\mathbf{F}}_{\alpha} - {\mathbf{A}}_{\alpha \beta} {\mathbf{u}}_{\beta} \right) +\underline{u} {\mathbf{A}}_{\alpha \alpha}^{-1}{\mathbf{A}}_{\alpha \beta} {\mathbf{1}}_{\beta}\\ \label{eq:matineq1} &= & {\mathbf{A}}_{\alpha \alpha}^{-1}{\mathbf{F}}_{\alpha} - {\mathbf{A}}_{\alpha \alpha}^{-1}{\mathbf{A}}_{\alpha \beta} \left({\mathbf{u}}_{\beta} -\underline{u}{\mathbf{1}}_{\beta}\right) \ge {\mathbf{0}}, \end{align}\tag{27}\] where we used 23 together with \({\mathbf{u}}_{\beta} -\underline{u}{\mathbf{1}}_{\beta}\ge {\mathbf{0}}\) and \({\mathbf{F}}_{\alpha} \ge {\mathbf{0}}\). Therefore, \[\min_D u_h \ge \underline{u} = \min_{\partial D}u_h^b.\] Moreover, if \({\mathbf{u}}_{\beta} -\underline{u}{\mathbf{1}}_{\beta} \gneqq {\mathbf{0}}\) or \({\mathbf{F}}_{\alpha}\gneqq {\mathbf{0}}\), then \(-{\mathbf{A}}_{\alpha \alpha}^{-1}{\mathbf{A}}_{\alpha \beta} <{\mathbf 0}\) together with \({\mathbf{A}}_{\alpha \alpha}^{-1}> {\mathbf{0}}\) and 27 imply that \({\mathbf{u}}_{\alpha}-\underline{u}{\mathbf{1}}_{\alpha} > {\mathbf{0}}\), showing that \(u_h\) cannot attain the global minimum \(\underline{u}\) in the interior of \(D\). So if \(u_h\) attains the value \(\underline{u}\) in the interior of \(D\), we have we must have \({\mathbf{u}}_{\beta} -\underline{u}{\mathbf{1}}_{\beta} = {\mathbf{0}}\) and \({\mathbf{F}}_{\alpha} = {\mathbf{0}}\). Now 27 implies \({\mathbf{u}}_{\alpha}-\underline{u}{\mathbf{1}}_{\alpha} = {\mathbf{0}}\), proving that \(u_h\) is constant. This shows that \({\mathcal{S}}^D_h\) satisfies sDMP-A.

To prove the reverse implication, assume \({\mathcal{S}}^D_h\) satisfies sDMP-A. For an internal node \(i\in \alpha\) define the source \(f^{(i)}\) so that \(\left< f^{(i)} , \varphi_j \right> = \delta_{ij}\), i.e., \(f^{(i)}\) is represented by the basis element \({\mathbf{e}}_i\); let \(u_h^b\equiv 0\). Due to sDMP-A, the vector \({\mathbf{u}}^{(i)}\) representing the solution \(u_h^{(i)} = {\mathcal{S}}^D_h(f^{(i)}, 0)\) satisfies \[\min {\mathbf{u}}^{(i)} \ge {\mathbf{0}}.\] In addition, if \({\mathbf{u}}^{(i)}\) has any zero entry, then \(u_h^{(i)}\) has a global minimum in the interior of \(D\), so by sDMP-A, \(u_h^{(i)}\equiv 0\), which, in turn, would imply \(f^{(i)}\equiv 0\). This shows \({\mathbf{u}}^{(i)} > {\mathbf{0}}\), meaning, the column corresponding to \(i\) in \({\mathbf{A}}_{\alpha \alpha}^{-1}\) is positive. Therefore, \({\mathbf{A}}_{\alpha \alpha}^{-1} > {\mathbf{0}}\). A similar argument, this time taking \(f\equiv 0\) and boundary vectors of the form \({\mathbf{u}}_{\beta} = {\mathbf{e}}_j\) with \(j\in \beta\) shows that \(-{\mathbf{A}}_{\alpha \alpha}^{-1}{\mathbf{A}}_{\alpha \beta} > {\mathbf{0}}\). ◻

Proof. Again, we restrict our attention to the case when \(f\) is nonnegative. The matrix form of 14 with \(\tilde{c} \ge 0\) is \[\label{eq:matforFEcpos} ({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha}){\mathbf{u}}_{\alpha} + ({\mathbf{A}}_{\alpha \beta} + {\mathbf{C}}_{\alpha \beta}){\mathbf{u}}_{\beta} = {\mathbf{F}}_{\alpha}.\tag{29}\] Using 25 , we have \[\begin{align} \label{eq:stiffCone} ({\mathbf{A}}_{\alpha \alpha} + {\mathbf{C}}_{\alpha \alpha}){\mathbf{1}}_{\alpha} + ({\mathbf{A}}_{\alpha \beta} + {\mathbf{C}}_{\alpha \beta}){\mathbf{1}}_{\beta} = {\mathbf{C}}_{\alpha \alpha} {\mathbf{1}}_{\alpha} + {\mathbf{C}}_{\alpha \beta}{\mathbf{1}}_{\beta}. \end{align}\tag{30}\] After multiplying 30 by a scalar \(r\) and subtracting from 29 we obtain \[\begin{align} ({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})({\mathbf{u}}_{\alpha} - r {\mathbf{1}}_{\alpha}) + ({\mathbf{A}}_{\alpha \beta} + {\mathbf{C}}_{\alpha \beta})({\mathbf{u}}_{\beta} - r {\mathbf{1}}_{\beta}) = {\mathbf{F}}_{\alpha} - r ({\mathbf{C}}_{\alpha \alpha} {\mathbf{1}}_{\alpha} + {\mathbf{C}}_{\alpha \beta}{\mathbf{1}}_{\beta}). \end{align}\] Hence, \[\begin{align} \label{eq:stiffConeB} {\mathbf{u}}_{\alpha} - r {\mathbf{1}}_{\alpha} &=& -\: ({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1}({\mathbf{A}}_{\alpha \beta} + {\mathbf{C}}_{\alpha \beta})({\mathbf{u}}_{\beta} - r {\mathbf{1}}_{\beta}) \\ \nonumber&&{+}\; ({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1}{\mathbf{F}}_{\alpha} \\ \nonumber&&-\;r({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1} ({\mathbf{C}}_{\alpha \alpha} {\mathbf{1}}_{\alpha} + {\mathbf{C}}_{\alpha \beta}{\mathbf{1}}_{\beta}). \end{align}\tag{31}\] First we assume the conditions 28 hold. Then \[({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1}{\mathbf{F}}_{\alpha}\ge 0.\] Let \(\underline{u} = \min u_h^b = \min {\mathbf{u}}_\beta\). If \(\underline{u} \ge 0\), then \({\mathbf{u}}_{\beta}^- = {\mathbf{0}}\). By taking \(r=0\) in 31 we get \[\begin{align} {\mathbf{u}}_{\alpha} \ge \overbrace{-\: ({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1}({\mathbf{A}}_{\alpha \beta} + {\mathbf{C}}_{\alpha \beta})}^{> {\mathbf{0}}} {\mathbf{u}}_{\beta} + \overbrace{({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1}}^{>{\mathbf 0}}{\mathbf{F}}_{\alpha}\ge {\mathbf{0}} , \end{align}\] which shows that \(\min {\mathbf{u}}_{\alpha} \ge 0 = -\max {\mathbf{u}}_{\beta}^-\). If one of the coordinates of \({\mathbf{u}}_{\alpha}\) is 0, then we must have \({\mathbf{u}}_{\beta}={\mathbf{0}}\) and \({\mathbf{F}}_{\alpha} = {\mathbf{0}}\); from 29 it follows that \({\mathbf{u}}_{\alpha} ={\mathbf{0}}\).

If \(\underline{u} < 0\), then we take \(r=\underline{u}\) in 31 , and we get \[\begin{align} {\mathbf{u}}_{\alpha} - \underline{u} {\mathbf{1}}_{\alpha} &=& \overbrace{-\: ({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1}({\mathbf{A}}_{\alpha \beta} + {\mathbf{C}}_{\alpha \beta})}^{> {\mathbf{0}}} \overbrace{({\mathbf{u}}_{\beta} - \underline{u} {\mathbf{1}}_{\beta})}^{\ge {\mathbf{0}}} \\ &&-\;\underline{u}\underbrace{({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1} ({\mathbf{C}}_{\alpha \alpha} {\mathbf{1}}_{\alpha} + {\mathbf{C}}_{\alpha \beta}{\mathbf{1}}_{\beta})}_{\ge {\mathbf{0}}} + \underbrace{({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1}}_{>{\mathbf 0}}{\mathbf{F}}_{\alpha}\ge {\mathbf{0}}, \end{align}\] where we used \(-\underline{u} > 0\) and \({\mathbf{C}}\ge {\mathbf{0}}\). Therefore, \({\mathbf{u}}_{\alpha} - \underline{u} {\mathbf{1}}_{\alpha}\ge {\mathbf{0}}\), showing that \[\min {\mathbf{u}}_{\alpha} \ge \underline{u} = -\max {\mathbf{u}}_{\beta}^-.\] If \({\mathbf{u}}_{\alpha} - \underline{u} {\mathbf{1}}_{\alpha}\) has a coordinate that is zero, then we have \[({\mathbf{u}}_{\beta} - \underline{u} {\mathbf{1}}_{\beta}) = {\mathbf{0}},\;\;{\mathbf{F}}_{\alpha} = {\mathbf{0}},\;\; \underline{u} ({\mathbf{C}}_{\alpha \alpha} {\mathbf{1}}_{\alpha} + {\mathbf{C}}_{\alpha \beta}{\mathbf{1}}_{\beta}) = {\mathbf{0}}.\] By applying 31 with \(r=\underline{u}\) we get \(({\mathbf{u}}_{\alpha} - \underline{u} {\mathbf{1}}_{\alpha}) = {\mathbf{0}}\), showing that \({\mathbf{u}}\) is constant. Hence, we proved sDMP-B holds.

The proof of the reverse implication is similar to that in Lemma 2. ◻

Proof. Assume \(f\) is nonnegative. Since the matrix form of 14 is given, as in Lemma 3, by 29 , the first condition in 28 ensures that the linear equation 29 has a unique solution; hence, the problem is well-posed. Let \(\underline{u} = \min u_h^b = \min {\mathbf{u}}_\beta\). We assume \(\underline{u} \ge 0\). By taking \(r=\underline{u}\) in 31 we get \[\begin{align} {\mathbf{u}}_{\alpha} - \underline{u} {\mathbf{1}}_{\alpha} &=& \overbrace{-\: ({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1}({\mathbf{A}}_{\alpha \beta} + {\mathbf{C}}_{\alpha \beta})}^{> {\mathbf{0}}} \overbrace{({\mathbf{u}}_{\beta} - \underline{u} {\mathbf{1}}_{\beta})}^{\ge {\mathbf{0}}} \\ &&-\;\underline{u}\underbrace{({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1} ({\mathbf{C}}_{\alpha \alpha} {\mathbf{1}}_{\alpha} + {\mathbf{C}}_{\alpha \beta}{\mathbf{1}}_{\beta})}_{\le {\mathbf{0}}} + \underbrace{({\mathbf{A}}_{\alpha \alpha} +{\mathbf{C}}_{\alpha \alpha})^{-1}}_{>{\mathbf 0}}{\mathbf{F}}_{\alpha}\ge {\mathbf{0}}, \end{align}\] where we used \(\underline{u} \ge 0\) and \({\mathbf{C}}\le {\mathbf{0}}\). Therefore, \({\mathbf{u}}_{\alpha} - \underline{u} {\mathbf{1}}_{\alpha}\ge {\mathbf{0}}\); thus, \(\min {\mathbf{u}}_{\alpha} \ge \underline{u}\). If \({\mathbf{u}}_{\alpha} - \underline{u} {\mathbf{1}}_{\alpha}\) has a coordinate that is zero, then we have \[({\mathbf{u}}_{\beta} - \underline{u} {\mathbf{1}}_{\beta}) = {\mathbf{0}},\;\;{\mathbf{F}}_{\alpha} = {\mathbf{0}},\;\; \underline{u} ({\mathbf{C}}_{\alpha \alpha} {\mathbf{1}}_{\alpha} + {\mathbf{C}}_{\alpha \beta}{\mathbf{1}}_{\beta}) = {\mathbf{0}}.\] By applying 31 with \(r=\underline{u}\) we get \(({\mathbf{u}}_{\alpha} - \underline{u} {\mathbf{1}}_{\alpha}) = {\mathbf{0}}\), showing that \({\mathbf{u}}\) is constant. Hence, we proved that \({\mathcal{S}}^D_h\) satisfies sDMP-R. ◻

Proof. We begin by defining a sequence of families of consistent discrete solution operators. This family would, in principle, be the finite element discretization of 1 2 with \({c}(x,u) = -\varepsilon u\) and \(\varepsilon>0\); however, neither the continuous, nor the discrete versions are guaranteed to be well-posed for a sufficiently rich family of \(\varepsilon>0\). In order to circumvent this problem, we restrict our attention to discrete solution operators, and only to a certain sequence \((\varepsilon_n)_{n\in \mathbb{N}}\) that converges to 0.

To define our solution operators first we note that the triangulation \({\mathcal{T}}_h\) is fixed. For every discrete subdomain \(E\subseteq D\) denote its sets of interior and boundary vertices by \(\alpha=\alpha(E)\) and \(\beta=\beta(E)\), respectively. Define the number \[\begin{align} \label{eq:setsigma} \delta = \min\bigcup_{E\subseteq D}\sigma({\mathbf{M}}_{\alpha(E) \alpha(E)}^{-1}{\mathbf{A}}_{\alpha(E) \alpha(E)}), \end{align}\tag{37}\] with the union being taken over all the discrete subdomains, and \(\sigma({\mathbf{B}})\) denoting the spectrum of a matrix \({\mathbf{B}}\). The expression in 37 is well defined because the set on the right-hand side above is finite. Moreover, since both \({\mathbf{M}}_{\alpha(E) \alpha(E)}\) and \({\mathbf{A}}_{\alpha(E) \alpha(E)}\) are symmetric positive definite matrices, the set in 37 lies in \((0,\infty)\), implying that \(\delta>0\). Let \(\varepsilon_n\) be a sequence so that \(0<\varepsilon_n<\delta\) and \(\lim_{n\to\infty }\varepsilon_n = 0\). Define \[\begin{align} {\mathbf{B}}^{(n)}_i = {\mathbf{A}}_{\alpha_i \alpha_i} - \varepsilon_n {\mathbf{M}}_{\alpha_i \alpha_i}. \end{align}\] Since \(\varepsilon_n\notin \sigma({\mathbf{M}}_{\alpha(E) \alpha(E)}^{-1}{\mathbf{A}}_{\alpha(E) \alpha(E)})\), it follows that \({\mathbf{B}}^{(n)}_i\) is invertible. Due to the continuity of matrix inversion we have the following: for all \(i\in I\) \[\begin{align} \lim_{n\to \infty} ({\mathbf{B}}^{(n)}_i)^{-1} = {\mathbf{A}}^{-1}_{\alpha_i \alpha_i} > {\mathbf{0}}. \end{align}\] Hence, there exists \(n_0\in \mathbb{N}\) so that for \[\label{eq:Bngzero} \forall n\ge n_0,\;\;\forall i\in I,\;\;({\mathbf{B}}^{(n)}_i)^{-1} > {\mathbf{0}}.\tag{38}\] It is assumed that \(\forall k\in \beta_i\), \(\exists j\in \alpha_i\) so that the boundary vertex \(P_k\) (of \(E_i\)) is connected to the interior vertex \(P_j\); this translates into \(({\mathbf{M}}_{\alpha_i \beta_i})_{jk}>0\), showing that every column of \({\mathbf{M}}_{\alpha_i \beta_i}\) has at least one positive entry. Since \({\mathbf{M}}_{\alpha_i \beta_i} \ge {\mathbf{0}}\), 38 implies that \[\begin{align} \label{eq:BnMnlzero1} \forall n\ge n_0,\;\;\forall i\in I,\;\; ({\mathbf{A}}_{\alpha_i \alpha_i} - \varepsilon_n {\mathbf{M}}_{\alpha_i \alpha_i})^{-1} {\mathbf{M}}_{\alpha_i \beta_i} > {\mathbf{0}}. \end{align}\tag{39}\] Therfore, it follows from 39 , 36 , and \(\varepsilon_n>0\) that \[\begin{align} \label{eq:BnMnlzero2} \;\;\;\forall n\ge n_0,\;\; \forall i\in I,\;\;({\mathbf{A}}_{\alpha_i \alpha_i} - \varepsilon_n {\mathbf{M}}_{\alpha_i \alpha_i})^{-1} ({\mathbf{A}}_{\alpha_i \beta_i}-\varepsilon_n{\mathbf{M}}_{\alpha_i \beta_i}) < {\mathbf{0}}. \end{align}\tag{40}\]

Now we construct the family of consistent solution operators, namely we consider the following sequence of discrete variational problems: given a discrete subdomain \(E\subseteq D\), find \(u_h\in V^h(\overline{E})\) of the form \(u_h = u_h^0+ u_h^b\) with \(u_h^0\in V_0^h(\overline{E})\), so that \[\begin{align} \label{eq:weakelldiscnegeps} a_E(u_h^0+u_h^b,v) -\varepsilon_n\left( u_h^0+u_h^b , v \right)_E=\left< f , v \right>_E,\;\;\forall v\in V_0^h(\overline{E}). \end{align}\tag{41}\] If \(\alpha=\alpha(E)\) and \(\beta=\beta(E)\), then 41 is formulated in matrix form as \[\begin{align} \label{eq:weakelldiscnegepsmat} ({\mathbf{A}}_{\alpha \alpha}-\varepsilon_n{\mathbf{M}}_{\alpha \alpha}){\mathbf{u}}_{\alpha} + ({\mathbf{A}}_{\alpha \beta} -\varepsilon_n{\mathbf{M}}_{\alpha \beta}){\mathbf{u}}_{\beta} = {\mathbf{F}}_{\alpha}. \end{align}\tag{42}\] The choice of \(\varepsilon_n\) ensures that 41 is well posed (uniquely solvable), hence this gives rise to a family of solution operators \({\mathcal{S}}^E_{h,n}: (V_0^h(\overline{E}))^*\times V^h(\partial E) \to V^h(E)\), so that \({\mathcal{S}}^E_{h,n}(f,u_{\beta}) = u = u_{\alpha}+u_{\beta}\), with \({\mathbf{u}}_{\alpha}\) and \({\mathbf{u}}_{\beta}\) satisfying 42 being the vector representations of \({u}_{\alpha}\) and \({u}_{\beta}\), respectively. Hence, Lemma 4 and 40 imply that for all \(i\in I\) and \(n\ge n_0\), \({\mathcal{S}}^{E_i}_{h,n}\) satisfies sDMP-R. Due to the variational formulation 41 , we can use an argument similar to that in Lemma 1 to show that for all \(n\ge n_0\), the family \(({\mathcal{S}}^E_{h,n})_{E\subseteq D}\) is consistent, as in Definition 1. Theorem 4(R) now applies to show that for all \(n\ge n_0\), \({\mathcal{S}}^{D}_{h,n}\) satisfies sDMP-R.

For the final step let \(f \in (V_0^h(\overline{D}))^*\) be nonnegative and \(u_h={\mathcal{S}}^{D}_{h}(f,u^b)\) be the solution of 14 with \(c\equiv 0\). Denote \[\begin{align} \underline{u} = \min_{\overline{D}} u_h,\;\;\; v_h = u_h-\underline{u}+1,\;\;\;\mathrm{and}\;\; v_h^b = u_h^b-\underline{u}+1. \end{align}\] Since \(c\equiv 0\), we have \(v_h={\mathcal{S}}^{D}_{h}(f,v_h^b)\), because the pair \({\mathbf{v}}_{\alpha}, {\mathbf{v}}_{\beta}\) given by \[{\mathbf{v}}_{\alpha} = {\mathbf{u}}_{\alpha}-(\underline{u}-1){\mathbf{1}}_{\alpha},\;\;\; {\mathbf{v}}_{\beta} = {\mathbf{u}}_{\beta}-(\underline{u}-1){\mathbf{1}}_{\beta},\] where \(\alpha=\alpha(D)\), \(\beta=\beta(D)\) satisfy 24 . This way we ensure \(v_h \ge 1\) (the number 1 is not special, all that matters is that \(1>0\)). Define \(v_{h,n}={\mathcal{S}}^{D}_{h,n}(f,v_h^b)\). Using 24 and 42 , we get \[\lim_{n\to\infty}v_{h,n} = v_h \ge 1\;\;\mathrm{in}\;\;V^h(D)\;\;\mathrm{(pointwise\;or\;any\;norm)},\] showing that \(\exists n_1\ge n_0\) so that for \(n\ge n_1\) we have \(v_{h,n}>0\). Since \({\mathcal{S}}^{D}_{h,n}\) satisfies sDMP-R, we get \[\label{eq:sdmprvn} \min_{\overline{D}} v_{h,n} = \min_{\partial{D}} v_{h,n} = \min_{\partial{D}} v^b_{h},\;\;\forall n\ge n_1.\tag{43}\] After passing to the limit in 43 we get \[\label{eq:sdmpru} \min_{\overline{D}} u_{h} = \min_{\overline{D}} v_{h} + \underbar{u}-1 = \min_{\partial{D}} v^b_{h} + \underbar{u}-1 = \min_{\partial{D}} u^b_{h},\tag{44}\] showing that \({\mathcal{S}}^{D}_{h}\) satisfies wDMP-A. ◻

Proof. Let \(E\subseteq D\) a connected discrete subdomain. For every vertex \(P_i\) interior to \(E\) we consider the set \(E_i = \bigcup \{T\in{\mathcal{T}}_h\; :\;P_i\in T\}\), which we call the star around \(P_i\). The boundary vertices of \(E_i\) are the adjacent vertices of \(P_i\), which is the only interior node of \(E_i\). In the terminology of Lemma 2, \(\alpha = \{i\}\) and \(\beta=\{j\;:\;{\mathbf{A}}_{ij}<0\}\), and the conditions 23 are satisfied in a trivial manner. Hence, \({\mathcal{S}}^{E_i}_h\) satisfies sDMP-A. Corollary 1 implies that \({\mathcal{S}}^{E}_h\) satisfies sDMP-A. ◻

Proof. For \(v\in V_0^h(\overline{D})\), the supports of \(v\) and \(\phi\) intersect in a set of measure zero, the third edge of \(T\). ◻

Proof. For (i), the hypothesis implies that each interior vertex \(P\) lies either in the interior of a patch \(F_P\) that is a rescaled version of \(G_{\ell}(\theta)\), or in a star \(E_P\) (the union of the triangles having \(P\) as vertex) that has only acute angles. In the former case, since the stiffness matrix in 2D is independent of the scale, Lemma 6 shows \({\mathbf{A}}^{-1}_{\ell} = {\mathbf{A}}^{-1}_{\ell}(\theta) > {\mathbf{0}}\), hence the same holds for the stiffness matrix \({\mathbf{A}}_P\) associated with the interior vertices of \(F_P\), because \({\mathbf{A}}_P\) is identical to \({\mathbf{A}}_{\ell}\). Now let \({\mathbf{A}}^b_{\ell}={\mathbf{A}}^b_{\ell}(\theta)\) be the stiffness matrix connecting interior and boundary entries in \(G_{\ell}(\theta)\), and \({\mathbf{A}}^b_P\) be its analogue on \(F_P\). Since the boundary layer of \(F_P\) contains only edges that satisfy 46 , we have \(({\mathbf{A}}^b_P)_{ij}<0\) for every edge \(\overline{P_i P_j}\) with \(P_i\in Int(F_P)\) and \(P_j\in \partial F_P\). This shows that \({\mathbf{A}}_P^{-1} {\mathbf{A}}^b_P < {\mathbf{0}}\), that is, 32 is verified for \(F_P\). For the latter case, when the patch is the star \(E_P\), the same argument as in Theorem 7 shows that 32 is verified on \(E_P\) as well. Corollary 2(i) shows that \({\mathcal{S}}^{\tilde{D}}_h\) satisfies sDMP-A. As in Example 1, adding the two corner triangles leads to \({\mathcal{S}}^{{D}}_h\) satisfying wDMP-A.
The argument for (ii) lies in the scaling of the mass matrix compared to that of the stiffness matrix. Assume a point \(P\) lies in the interior of a patch \(F_P\) that is a rescaled version of \(G_{\ell}(\theta)\), as above. Denote by \({\mathbf{M}}_{\ell}={\mathbf{M}}_{\ell}(\theta)\) the mass matrix associated with the interior vertices of the reference mesh \(G_{\ell}(\theta)\). Due to scaling (see 49 ), the associated mass matrix for \(F_P\) is \(K h^2{\mathbf{M}}_{\ell}\), with \(K\) a dimensionless constant depending on \(\theta\) and \(\ell\). Similarly, if \({\mathbf{M}}^b_{\ell}={\mathbf{M}}^b_{\ell}(\theta)\) is the mass matrix connecting interior and boundary entries on \(G_{\ell}(\theta)\), then the analogous matrix on \(F_P\) is \(K h^2{\mathbf{M}}^b_{\ell}\). Hence, the interior and boundary reaction matrices \({\mathbf{C}}_P\) and \({\mathbf{C}}^b_P\), respectively, associated with \(F_P\) (see 22 ) satisfy \[\label{eq:masstozero} {\mathbf{0}}\le {\mathbf{C}}_P \le K \tilde{c}_{\max} h^2{\mathbf{M}}_{\ell},\;\;\;{\mathbf{0}}\le {\mathbf{C}}^b_P \le K \tilde{c}_{\max} h^2{\mathbf{M}}^b_{\ell}.\tag{53}\] Consequently \[\begin{align} \lim_{h\to 0} ({\mathbf{A}}_P + {\mathbf{C}}_P)^{-1} = {\mathbf{A}}_P^{-1} > {\mathbf{0}}, \end{align}\] and\[\begin{align} \lim_{h\to 0} ({\mathbf{A}}_P + {\mathbf{C}}_P)^{-1}({\mathbf{A}}^b_P + {\mathbf{C}}^b_P) = {\mathbf{A}}_P^{-1}{\mathbf{A}}^b_P < {\mathbf{0}}. \end{align}\] Hence, there exists \(h_{\max}>0\) depending on \(\tilde{c}_{\max}\) and \(K\) (which depends on \(\theta\) and \(k\) - the largest index of the defects), so that for \(0<h<h_{\max}\) we have \[\begin{align} \label{eq:Cconditionsex1} ({\mathbf{A}}_P + {\mathbf{C}}_P)^{-1} > {\mathbf{0}},\;\;\mathrm{and}\;\; ({\mathbf{A}}_P + {\mathbf{C}}_P)^{-1}({\mathbf{A}}^b_P + {\mathbf{C}}^b_P) < {\mathbf{0}}. \end{align}\tag{54}\] If, instead, the star \(E_P\) around \(P\) has only acute angles, similar inequalities hold if we replace \(F_P\) with \(E_P\). By Corollary 2(ii), satisfies \({\mathcal{S}}^{\tilde{D}}_h\) satisfies sDMP-B. The argument shown in Example 1 can be replicated to show that when adding back the two triangles at the horizontal extremes of \(\tilde{D}\), the global solution operator \({\mathcal{S}}^{{D}}_h\) satisfies wDMP-B. ◻

Proof. The statement follows from Corollary 3 and Theorem 6. The hypotheses ensure that every vertex \(V\) lies in the interior of a patch similar to \(E\) in the top left corner of Fig. 12, or the star around \(V\) is similar to a standard star in Example 1; in either case, the inequalities 36 are satisfied. ◻

Proof. The argument follows the same path as in Theorem 7, namely we prove that \({\mathcal{S}}^E_h\) satisfies the sDMP-B for any star \(E_P\) around a vertex \(P\). The conclusion will follow from Corollary 1.

Thus, it is sufficient to prove the prove sDMP-B assuming the mesh has only one interior point \(P_1\) and boundary points \(P_2, \dots P_N\), i.e., \(D=E_{P_1}\). In this case \(h \le \mathrm{diam}(D) \le 2 h\); note that the mesh size plays a role in the estimation. Let \(u_h={\mathcal{S}}_h^D(f,u_h^b)\), and denote \({\mathbf{F}}_i = f(P_i)\), \({\mathbf{U}}_i = u_h(P_i)\), for \(1\le i\le N\). In matrix terms, 15 is represented by the nonlinear scalar equation: find \({\mathbf{U}}_1 \in \mathbb{R}\) so that \[\begin{align} \label{eq:sldiscmat} {\mathbf{A}}_{11} {\mathbf{U}}_1 + \sum_{i=2}^N {\mathbf{A}}_{i1}{\mathbf{U}}_i + {\mathbf{M}}_{11} \:c(P_1, {\mathbf{U}}_1) + \sum_{i=2}^N {\mathbf{M}}_{i1}\:c(P_i, {\mathbf{U}}_i) = \sum_{i=1}^N {\mathbf{M}}_{i1} {\mathbf{F}}_i. \end{align}\tag{57}\] Since \(f\ge 0\), we have \[\begin{align} \label{eq:sldiscmatineq} {\mathbf{A}}_{11} {\mathbf{U}}_1 + \sum_{i=2}^N {\mathbf{A}}_{i1}{\mathbf{U}}_i + {\mathbf{M}}_{11} \:c(P_1, {\mathbf{U}}_1) + \sum_{i=2}^N {\mathbf{M}}_{i1}\:c(P_i, {\mathbf{U}}_i) \ge 0. \end{align}\tag{58}\] Note that \({\mathbf{A}}_{11}>0\), and \({\mathbf{M}}_{i1}>0\) for \(1\le i \le N\). We consider two cases:
Case 1: First assume that \({\mathbf{U}}_i\ge 0\) for \(2\le i\le N\), which implies that \[\label{eq:posClip} 0\le c(P_i, {\mathbf{U}}_i)\le L_c\: {\mathbf{U}}_i,\;\;2\le i\le N.\tag{59}\] In this case we have \[-\max_{\partial D}{u^-_h} = -\max_{2\le i\le N} {\mathbf{U}}_i^- = 0.\] If we had \({\mathbf{U}}_1<0\), then \(c(P_1, {\mathbf{U}}_1)\le 0\). Hence, \[\begin{align} \nonumber &&0<{\mathbf A}_{11} (-{\mathbf{U}}_1) - {\mathbf{M}}_{11} \:c(P_1, {\mathbf{U}}_1) \stackrel{\eqref{eq:sldiscmatineq}}{\le} \sum_{i=2}^N {\mathbf{A}}_{i1}{\mathbf{U}}_i + \sum_{i=2}^N {\mathbf{M}}_{i1}\:c(P_i, {\mathbf{U}}_i)\\ \nonumber &&\stackrel{\eqref{eq:posClip}}{\le} \sum_{i=2}^N ({\mathbf{A}}_{i1} + L_c{\mathbf{M}}_{i1}){\mathbf{U}}_i =\sum_{i=2}^N (-{\mathbf{A}}_{i1})(-1 + L_c{\mathbf{M}}_{i1}/(-{\mathbf{A}}_{i1})){\mathbf{U}}_i\\ \label{eq:case1ineq} &&\stackrel{\eqref{eq:stiffmatrixlimitsrat}}{\le} \sum_{i=2}^N (-{\mathbf{A}}_{i1})(-1 + L_c C_A h^2){\mathbf{U}}_i. \end{align}\tag{60}\] This is not possible if \(h<h_0 = \sqrt{(L_c C_A)^{-1}}\), for in that case we have \[\label{eq:case1ineq2} (-{\mathbf{A}}_{i1})(-1 + L_c C_A h^2) < 0,\;\;\forall i=2,\dots,N.\tag{61}\] Together with \({\mathbf{U}}_i\ge 0\) for \(2\le i\le N\), this renders the sum in 60 to be nonpositive, contradicting our assumption on \({\mathbf{U}}_1\). It remains that \({\mathbf{U}}_1\ge 0\). For the strong part, assume \({\mathbf{U}}_1 = 0\); then we also have \(c(P_1, {\mathbf{U}}_1) = 0\). Using the same line of arguments as above, it follows that the sum from 60 is nonnegative. However, since all the terms in 60 are nonpositive, they all must be zero, forcing \({\mathbf{U}}_i = 0\) for \(i=2,\dots,N\).
Case 2: We assume that at least one of the values \({\mathbf{U}}_i\), \(i=2, \dots N\) is negative. We partition \(I_N = \{2,\dots,N\} = {\mathcal{N}}_{-}\cup {\mathcal{N}}_{+}\), with \[\begin{align} {\mathcal{N}}_{-} = \{i\in I_N :\;{\mathbf{U}}_i< 0\},\;\;{\mathcal{N}}_{+} = \{i\in I_N :\;{\mathbf{U}}_i\ge 0\}. \end{align}\] This implies \(c(P_i, {\mathbf{U}}_i)\le 0\) for \(i\in {\mathcal{N}}_{-}\), and \(c(P_i, {\mathbf{U}}_i)\ge 0\) for \(i\in {\mathcal{N}}_{+}\). Note that \[\begin{align} \label{eq:acoeffineq} -{\mathbf{A}}_{11} - \sum_{i\in {\mathcal{N}}_{-}} {\mathbf{A}}_{i1} = \sum_{i\in {\mathcal{N}}_{+}} {\mathbf{A}}_{i1} \le 0, \end{align}\tag{62}\] because \(\sum_{i = 1}^N{\mathbf{A}}_{1i}=0\). Note that all the arguments still hold if \({\mathcal{N}}_+=\varnothing\). If \({\mathbf{U}}_1 < \min_{i\in {\mathcal{N}}_{-}} {\mathbf{U}}_i<0\), then cf. 58 we have \[\begin{align} & \sum_{i\in {\mathcal{N}}_{+}} {\mathbf{A}}_{i1}{\mathbf{U}}_i + {\mathbf{M}}_{11} \:c(P_1, {\mathbf{U}}_1) + \sum_{i=2}^N {\mathbf{M}}_{i1}\:c(P_i, {\mathbf{U}}_i) \ge -{\mathbf{A}}_{11} {\mathbf{U}}_1 + \sum_{i\in {\mathcal{N}}_{-}} (-{\mathbf{A}}_{i1}){\mathbf{U}}_i\\ & = \underbrace{(-{\mathbf{A}}_{11} -\sum_{i\in {\mathcal{N}}_{-}} {\mathbf{A}}_{i1}) {\mathbf{U}}_1}_{\ge 0} + \sum_{i\in {\mathcal{N}}_{-}} (-{\mathbf{A}}_{i1})({\mathbf{U}}_i-{\mathbf{U}}_1)\ge \sum_{i\in {\mathcal{N}}_{-}} (-{\mathbf{A}}_{i1})({\mathbf{U}}_i-{\mathbf{U}}_1) >0. \end{align}\] Therefore, \[\begin{align} \nonumber &0< {\mathbf{M}}_{11} \:c(P_1, {\mathbf{U}}_1) + \sum_{i\in {\mathcal{N}}_{+}} \left({\mathbf{A}}_{i1}{\mathbf{U}}_i+{\mathbf{M}}_{i1}\:c(P_i, {\mathbf{U}}_i) \right) + \sum_{i\in {\mathcal{N}}_{-}} {\mathbf{M}}_{i1}\:c(P_i, {\mathbf{U}}_i)\\ \label{eq:sldiscmatineqcase3} &\stackrel{\eqref{eq:stiffmatrixlimitsrat}}{\le} {\mathbf{M}}_{11} \:c(P_1, {\mathbf{U}}_1) + \sum_{i\in {\mathcal{N}}_{+}} (-{\mathbf{A}}_{i1})(-1 + L_c C_A h^2){\mathbf{U}}_i + \sum_{i\in {\mathcal{N}}_{-}} {\mathbf{M}}_{i1}\:c(P_i, {\mathbf{U}}_i). \end{align}\tag{63}\] Following 61 , all the terms in the sum of 63 are nonpositive for \(h<h_0\), thus contradicting the assumption on \({\mathbf{U}}_1\). Hence, \[\label{eq:minu0} {\mathbf{U}}_1 \ge \min_{i\in {\mathcal{N}}_{-}} {\mathbf{U}}_i = -\max_{\partial E_{P_1}} u_h^-.\tag{64}\] If equality holds in 64 , then we have \({\mathbf{U}}_i-{\mathbf{U}}_1 \ge 0\), \(c(P_i, {\mathbf{U}}_i)-c(P_i, {\mathbf{U}}_1) \ge 0\), and \(c(P_i, {\mathbf{U}}_1)\le 0\) for \(1\le i \le N\). Consequently, \[\begin{align} \label{eq:sldiscmatineq4} &0\le {\mathbf{A}}_{11}{\mathbf{U}}_1 + \sum_{i=2}^N {\mathbf{A}}_{i1}{\mathbf{U}}_i + \sum_{i=1}^N{\mathbf{M}}_{i1}\:c(P_i,{\mathbf{U}}_i) \\ &= \sum_{i=2}^N {\mathbf{A}}_{i1}({\mathbf{U}}_i-{\mathbf{U}}_1) + \sum_{i=1}^N{\mathbf{M}}_{i1}\:c(P_i,{\mathbf{U}}_1) + \sum_{i=2}^N{\mathbf{M}}_{i1}\:\left(c(P_i,{\mathbf{U}}_i) - c(P_i, {\mathbf{U}}_1)\right)\\ &\le \sum_{i=2}^N ({\mathbf{A}}_{i1}+ L_c {\mathbf{M}}_{i1})({\mathbf{U}}_i-{\mathbf{U}}_1) + \sum_{i=1}^N{\mathbf{M}}_{i1}\:c(P_i,{\mathbf{U}}_1)\\ &\le \sum_{i=2}^N (-{\mathbf{A}}_{i1})(-1 + L_c C_A h^2 )({\mathbf{U}}_i-{\mathbf{U}}_1) + \sum_{i=1}^N{\mathbf{M}}_{i1}\:c(P_i,{\mathbf{U}}_1). \end{align}\tag{65}\] Since all the terms in the sum above are nonpositive for \(h<h_0\), it follows that they are all zero, showing that \({\mathbf{U}}\) is constant and that \(c(P_i,{\mathbf{U}}_1) = 0\) for \(1\le i\le N\). ◻

Proof. Referring to Figures 10 and 14, note that \[\begin{align} \widetilde{\varphi}_B& =& {\varphi}_B + \widetilde{\varphi}_B(M) {\varphi}_M + \widetilde{\varphi}_B(N) {\varphi}_N + \widetilde{\varphi}_B(P) {\varphi}_P\\ & =& {\varphi}_B + \frac{1}{2} {\varphi}_M + \left(\frac{3}{4}+O(\alpha)\right) {\varphi}_N + \left(\frac{1}{4}+O(\alpha)\right) {\varphi}_P, \end{align}\] where we used the fact that \(N=N(\alpha) \to (1/4,0)\) as \(\alpha\to 0\), showing that \(\lim_{\alpha\to 0} \widetilde{\varphi}_B(N) = 3/4\). Due to the smoothness of the function \(\alpha\mapsto \widetilde{\varphi}_B(N(\alpha))\), we conclude that \(\widetilde{\varphi}_B(N) = 3/4 + O(\alpha)\). Similarly, \(\widetilde{\varphi}_B(P) = 1/4 + O(\alpha)\). Following the same line of arguments we have \[\begin{align} \widetilde{\varphi}_C& =& {\varphi}_C + \widetilde{\varphi}_C(M) {\varphi}_M + \widetilde{\varphi}_C(N) {\varphi}_N + \widetilde{\varphi}_C(P) {\varphi}_P\\ & =& {\varphi}_C + \frac{1}{2} {\varphi}_M + \left(\frac{1}{4}+O(\alpha)\right) {\varphi}_N + \left(\frac{3}{4}+O(\alpha)\right) {\varphi}_P, \end{align}\] and \[\begin{align} \widetilde{\varphi}_A& =& {\varphi}_A + \widetilde{\varphi}_A(M) {\varphi}_M + \widetilde{\varphi}_A(N) {\varphi}_N + \widetilde{\varphi}_A(P) {\varphi}_P\\ & =& {\varphi}_A + O(\alpha) {\varphi}_N + O(\alpha) {\varphi}_P. \end{align}\] If we redenote the basis functions (for the purpose of temporarily replacing the letter indices with numbers) \({{\mathcal{B}}} = \{{\psi}_i\}_{i=1,\dots,6}\) and \(\widetilde{{\mathcal{B}}} = \{\widetilde{\psi}_i\}_{i=1,\dots,6}\), then \[\begin{align} \psi_i = \sum_{j=1}^6 {\mathbf{E}}_{ij} \widetilde{\psi}_j,\;\;i=1,\dots,6,\;\;\mathrm{with} \end{align}\] \[\begin{align} {\mathbf{E}} = \begin{bmatrix} {}1{}& {}0{}& {}0{}& -\varphi_B(M)& -\varphi_B(N)& -\varphi_B(P)\\ 0& 1& 0& -\varphi_C(M)& -\varphi_C(N)& -\varphi_C(P)\\ 0& 0& 1& -\varphi_A(M)& -\varphi_A(N)& -\varphi_A(P)\\ 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1 \end{bmatrix} \end{align}\]

Therefore, \[\begin{align} {\mathbf{S}}_{ij} & =& a_E(\psi_j,\psi_i) = a_E\left(\sum_{k=1}^6{\mathbf{E}}_{jk}\widetilde{\psi}_k,\sum_{k=1}^6{\mathbf{E}}_{il}\widetilde{\psi}_l\right) =\sum_{k,l=1}^6{\mathbf{E}}_{jk} {\mathbf{E}}_{il}\: a_E\left(\widetilde{\psi}_k,\widetilde{\psi}_l\right)\\ & = & \sum_{k,l=1}^6{\mathbf{E}}_{jk} {\mathbf{E}}_{il} \widetilde{{\mathbf{S}}}_{lk} = \left({\mathbf{E}} \widetilde{{\mathbf{S}}}{\mathbf{E}}^T\right)_{ij}, \end{align}\] showing 70 . ◻

Proof. Using 70 we get \[\label{eq:smallstiffnessmatinv} {\mathbf{S}}^{-1} = {\mathbf{E}}^{-T}\widetilde{{\mathbf{S}}}^{-1}{\mathbf{E}}^{-1}.\tag{72}\] For computing \(\widetilde{{\mathbf{S}}}^{-1}\) we use the block inversion formula \[\label{eq:blockinverse} \widetilde{{\mathbf{S}}}^{-1} = \begin{bmatrix} {\mathbf{I}} & {\mathbf{0}} \\ {\mathbf{D}}^T & {\mathbf{I}} \end{bmatrix}\; \begin{bmatrix} {\widehat{{\mathbf{S}}}}^{-1} & {\mathbf{0}} \\ {\mathbf{0}} & {{\mathbf{C}}}^{-1} \end{bmatrix}\; \begin{bmatrix} {\mathbf{I}} & {\mathbf{D}} \\ {\mathbf{0}} & {\mathbf{I}} \end{bmatrix}\tag{73}\] with \(\widehat{{\mathbf{S}}} = {{\mathbf{A}}} - {{\mathbf{B}}} {{\mathbf{C}}}^{-1}{{\mathbf{B}}}^T\) being the Schur complement of \({{\mathbf{C}}}\), and \({\mathbf{D}} = -{{\mathbf{B}}}{{\mathbf{C}}}^{-1}\). Note that \(\alpha {\mathbf{C}} \approx {\mathbf{C}_0}\), hence \(\alpha^{-1}{\mathbf{C}}^{-1} \approx {\mathbf{C}_0}^{-1}\), showing that \[\begin{align} {\mathbf{C}}^{-1} = \alpha\left({\mathbf{C}_0}^{-1}+{\mathbf{O}}(\alpha)\right) = \alpha \left( \begin{bmatrix} 6 & 4 & 4\\ 4 & 7/2 & 5/2\\ 4 & 5/2 & 7/2 \end{bmatrix} +{\mathbf{O}}(\alpha)\right) . \end{align}\] Since \({\mathbf{B}} = {\mathbf{O}}(1)\) and \({\mathbf{C}}^{-1} = {\mathbf{O}}(\alpha)\), it follows that \[\begin{align} \widehat{{\mathbf{S}}} = {{\mathbf{A}}} - {{\mathbf{B}}} {{\mathbf{C}}}^{-1}{{\mathbf{B}}}^T = {{\mathbf{A}}} + {\mathbf{O}}(\alpha), \end{align}\] showing that \(\widehat{{\mathbf{S}}} \approx {{\mathbf{A}}}\). Hence, \[\begin{align} \label{eq:Ainv} \widehat{{\mathbf{S}}}^{-1} \approx {\mathbf{A}}^{-1} = \begin{bmatrix} 2/7 & 1/14 & 1/14\\ 1/14 & 15/56 & 1/56\\ 1/14 & 1/56 & 15/56 \end{bmatrix} > {\mathbf{0}}. \end{align}\tag{74}\] Therefore, using 73 we get \[\begin{align} \label{eq:blockinverseeplicit} &&\widetilde{{\mathbf{S}}}^{-1} = \begin{bmatrix} {\mathbf{I}} & {\mathbf{0}} \\ {\mathbf{D}}^T & {\mathbf{I}} \end{bmatrix}\; \begin{bmatrix} {\widehat{{\mathbf{S}}}}^{-1} & {\mathbf{0}} \\ {\mathbf{0}} & {{\mathbf{C}}}^{-1} \end{bmatrix}\; \begin{bmatrix} {\mathbf{I}} & {\mathbf{D}} \\ {\mathbf{0}} & {\mathbf{I}} \end{bmatrix} = \begin{bmatrix} {\widehat{{\mathbf{S}}}}^{-1} & {\widehat{{\mathbf{S}}}}^{-1}{\mathbf{D}} \\ {\mathbf{D}}^T {\widehat{{\mathbf{S}}}}^{-1} & {{\mathbf{C}}}^{-1} + {\mathbf{D}}^T {\widehat{{\mathbf{S}}}}^{-1}{\mathbf{D}} \end{bmatrix} \end{align}\tag{75}\] Note that \[\begin{align} {\mathbf{D}} &=& -{{\mathbf{B}}}{{\mathbf{C}}}^{-1} = -\alpha\begin{bmatrix} {1}/{2}& 0& 0\\{1}/{2}& 0& 0\\-{1}/{2}& 0& 0 \end{bmatrix} \left(\begin{bmatrix} 6 & 4 & 4\\ 4 & 7/2 & 5/2\\ 4 & 5/2 & 7/2 \end{bmatrix} +{\mathbf{O}}(\alpha)\right)\\ &=& \alpha \left(\begin{bmatrix} -3& -2& -2\\-3& -2& -2\\3& 2& 2 \end{bmatrix}+{\mathbf{O}}(\alpha)\right) = \alpha \left({\mathbf{D}}_0 + {\mathbf{O}}(\alpha)\right). \end{align}\] Therefore \[\begin{align} {{\mathbf{C}}}^{-1} + {\mathbf{D}}^T {\widehat{{\mathbf{S}}}}^{-1}{\mathbf{D}} = {\alpha}({\mathbf{C}}^{-1}_0+{\mathbf{O}}(\alpha)) + {\mathbf{O}}(\alpha^2) ={\alpha}({\mathbf{C}}^{-1}_0+{\mathbf{O}}(\alpha)). \end{align}\] Hence \[\begin{align} \label{eq:blockinverseeplicit2} &&\widetilde{{\mathbf{S}}}^{-1} = \begin{bmatrix} {{{\mathbf{A}}}}^{-1} + {\mathbf{O}}(\alpha) & \alpha {{{\mathbf{A}}}}^{-1}\left({\mathbf{D}}_0 + {\mathbf{O}}(\alpha)\right) \\ \alpha \left({\mathbf{D}}^T_0 + {\mathbf{O}}(\alpha)\right) {{{\mathbf{A}}}}^{-1} & {\alpha}({\mathbf{C}}^{-1}_0+{\mathbf{O}}(\alpha)) \end{bmatrix} \end{align}\tag{76}\] Putting together 71 , 72 , and 76 we get \[\begin{align} {\mathbf{S}}^{-1} &=& \begin{bmatrix} {{{\mathbf{A}}}}^{-1} & {{{\mathbf{A}}}}^{-1} {\mathbf{R}}_0\\ {\mathbf{R}}^T_0 {{\mathbf{A}}}^{-1}& {\mathbf{R}}^T_0 {{\mathbf{A}}}^{-1} {\mathbf{R}}_0 \end{bmatrix} +{\mathbf{O}}(\alpha) = {\mathbf{T}}_0 + {\mathbf{O}}(\alpha). \end{align}\] A direct computation shows that \[\begin{align} {\mathbf{A}}^{-1} {\mathbf{R}}_0 = \begin{bmatrix} 5/28 & 13/56 & 1/8\\ 19/112 & 27/224 & 7/32\\ 5/112& 13/224 & 1/32 \end{bmatrix} \end{align}\] and \[\begin{align} {\mathbf{R}}^T_0{\mathbf{A}}^{-1} {\mathbf{R}}_0 = \begin{bmatrix} 39/224 & 79/448 & 11/64\\ 79/448 & 183/896 & 19/128\\ 11/64 & 19/128 & 25/128 \end{bmatrix}, \end{align}\] showing, together with 74 , that \({\mathbf{T}}_0 >{\mathbf 0}\). Since \({\mathbf{S}}^{-1} = {\mathbf{T}}_0 +{\mathbf{O}}(\alpha)\), we have \[\lim_{\alpha\to 0} {\mathbf{S}}^{-1} = {\mathbf{T}}_0.\] Note that \({\mathbf{T}}_0\) is singular, and has rank 3, hence is rank-3 deficient. The computation has been validated numerically against stiffness matrix computations using standard finite element codes. ◻

Proof. The fact that the \({{\mathbf{A}}}\) block has the form 68 is a simple consequence of the formulas 45 and 47 , and are well known for a uniform grid. Due to the fact \(\widetilde{\varphi}_A\), \(\widetilde{\varphi}_B\), \(\widetilde{\varphi}_C\) are linear on \(\Delta ABC\), and \(\varphi_N, \varphi_P\) vanish on \(\partial (\Delta ABC)\), we have \[\begin{align} a_E(\widetilde{\varphi}_X, \varphi_Y) = 0,\;\;\forall X \in \{A, B, C\},\;\; Y \in \{N, P\}. \end{align}\] (see also Lemma 5.5 in [17]), thus justifying all the six zero-entries in the matrix \({{\mathbf{B}}}\). The remaining non-trivial entries are associated with the set \(\{\varphi_M, \varphi_N, \varphi_P\}\).

Entries related to \({\varphi}_N, {\varphi}_P\): This case has been discussed in Lemma 5.6 in [17]; for completeness we review the computation, which is simplified because we are interested primarily in the case when \(0< \alpha \ll 1\). Refer to Figure 14 for notation. We have \[\begin{align} a_E(\varphi_N, \varphi_N) &= &\int_{{\mathcal{T}}_1\cup {\mathcal{T}}_2\cup {\mathcal{T}}_6\cup{\mathcal{T}}_7} |\nabla\varphi_N|^2 \\ &= &\frac{\sin(\pi-2 \alpha)}{2 \sin^2\alpha} + \frac{\sin\alpha}{2 \sin(\pi-2 \alpha) \sin\alpha} + \int_{{\mathcal{T}}_6\cup{\mathcal{T}}_7} |\nabla \varphi_N|^2\\ & = & \frac{\cos \alpha}{\sin\alpha} + \frac{1}{2 \sin(2 \alpha)} + O(1) = \frac{1}{\alpha}\left(\frac{5}{4} + O(\alpha) \right). \end{align}\] Similarly, \[\begin{align} a_E(\varphi_P, \varphi_P) &= & \frac{1}{\alpha}\left(\frac{5}{4} + O(\alpha) \right). \end{align}\] Furthermore, if \(\gamma\) denotes the angle \(\sphericalangle NAP\), then \[\begin{align} a_E(\varphi_N, \varphi_P)& = & -\frac{\sin(\gamma+\pi-2\alpha)}{2 \sin \gamma \sin(\pi-2\alpha)} = -\frac{\sin(2\alpha-\gamma)}{2 \sin \gamma \sin(2\alpha)} \\ &= & -\frac{\cos \gamma}{2 \sin \gamma} +\frac{\cos(2\alpha)}{2 \sin(2\alpha)} = \frac{1}{\alpha}\left(\frac{1}{4} + O(\alpha) \right) \end{align}\] because \(\lim_{\alpha\to 0} \gamma(\alpha) = \gamma_0\in (0,\pi/4)\), showing that the expression involving \(\gamma\) is bounded as \(\alpha\to 0\).

Entries related to \({\varphi}_M\): The quantities \(a_E(\varphi_M, \widetilde{\varphi}_X)\) with \(X \in \{A, B, C\}\) cannot be computed using 45 , because the latter are not nodal basis functions with respect to the finer mesh. First note that \(\nabla \widetilde{\varphi}_A = e_2\), showing that \[\begin{align} a_E(\varphi_M, \widetilde{\varphi}_A)& = & \int_{{\mathcal{T}}_1\cup {\mathcal{T}}_2\cup {\mathcal{T}}_3} \partial_y\varphi_M. \end{align}\] Let \(g:[0,1]\to E\) be the piecewise linear function whose graph is represented by the contour \(BNPC\). Using Fubini’s theorem and the fact that \(\varphi_M(g(x)) = 0\), we obtain \[\begin{align} \int_{{\mathcal{T}}_1\cup {\mathcal{T}}_2\cup {\mathcal{T}}_3} \partial_y\varphi_M & = & \int_0^1 dx \int_0^{g(x)}\partial_y\varphi_M(x,y) = - \int_0^1 \varphi_M(x,0) dx = -\frac{1}{2}, \end{align}\] because \(\varphi_M(x,0)\) is a one-dimensional, piecewise linear nodal basis function on \(\overline{BC}\). Therefore \(\widetilde{{S}}_{34} = \widetilde{{S}}_{43}=-1/2\). For similar reasons, if \(A'=(1,-1)\) denotes the vertex on \(\partial E\) directly below \(C\) (see Figure 14), and \(\widetilde{\varphi}_{A'}\) is the coarse nodal basis function associated with \(A'\), then \[\begin{align} a_E(\varphi_M, \widetilde{\varphi}_{A'}) &=& \int_{{\mathcal{T}}_4\cup {\mathcal{T}}_5} \nabla {\varphi}_{M} \cdot \nabla \widetilde{\varphi}_{A'} = -\frac{1}{2}. \end{align}\]

Since \(\mathrm{supp}(\varphi_M) = \bigcup_{i=1}^5 {\mathcal{T}}_i \subseteq\mathrm{supp}(\widetilde{\varphi}_B)\cap \mathrm{supp}(\widetilde{\varphi}_C)\), we have \[\begin{align} a_E(\varphi_M, \widetilde{\varphi}_X)& = & \int_{{\mathcal{T}}_1\cup\dots\cup {\mathcal{T}}_5} \nabla \varphi_M\cdot \nabla \widetilde{\varphi}_X,\;\;X\in \{B,C\}. \end{align}\] Note that on \(D_M^{(1)} = {\mathcal{T}}_1\cup {\mathcal{T}}_2 \cup {\mathcal{T}}_3\) we have \(\nabla \widetilde{\varphi}_C = e_1\), while on \(D_M^{(2)} = {\mathcal{T}}_4\cup {\mathcal{T}}_5\) we have \(\nabla \widetilde{\varphi}_B = -e_1\). Hence, after using Fubini’s theorem \[\begin{align} \label{eq:Dm1} &&\int_{D_M^{(1)}} \nabla \varphi_M \cdot \nabla \widetilde{\varphi}_C = \int_{D_M^{(1)}} \partial_x \varphi_M = \int_0^{\frac{1}{4}\tan \alpha } dy \int_{D_{M,y}^{(1)}} \partial_x \varphi_M \;dx = 0, \end{align}\tag{77}\] where \(D_{M,y}^{(1)}\) is the horizontal section of \(D_{M}^{(1)}\) at level \(y\), and we use the fact that \(\varphi_M\) is zero at the end points of each horizontal section of \(D_{M}^{(1)}\). Similarly, \[\begin{align} \label{eq:Dm2} && \int_{D_M^{(2)}} \nabla \varphi_M \cdot \nabla \widetilde{\varphi}_B = -\int_{D_M^{(2)}} \partial_x \varphi_M = -\int_{-1}^{0} dy \int_{D_{M,y}^{(2)}} \partial_x \varphi_M \; dx = 0, \end{align}\tag{78}\] for the same reason. It remains that \[\begin{align} \label{eq:MBC} && a_E(\varphi_M, \widetilde{\varphi}_B) = \int_{D_M^{(1)}} \nabla \varphi_M\cdot \nabla \widetilde{\varphi}_B,\;\;\; a_E(\varphi_M, \widetilde{\varphi}_C) = \int_{D_M^{(2)}} \nabla \varphi_M\cdot \nabla \widetilde{\varphi}_C. \end{align}\tag{79}\]

On \(\Delta ABC\) we have \(\widetilde{\varphi}_B + \widetilde{\varphi}_C + \widetilde{\varphi}_A \equiv 1\), thus (implicitly) this holds on \(D_M^{(1)}\). So \[\begin{align} 0 &= & \int_{D_M^{(1)}} \nabla \varphi_M \cdot \nabla (\widetilde{\varphi}_B + \widetilde{\varphi}_C + \widetilde{\varphi}_A) \stackrel{\eqref{eq:Dm1}}{=} \int_{D_M^{(1)}} \nabla \varphi_M \cdot \nabla \widetilde{\varphi}_B + \int_{D_M^{(1)}} \nabla \varphi_M \cdot \nabla \widetilde{\varphi}_A\\ & \stackrel{\eqref{eq:MBC}}{=} & a_E(\varphi_M, \widetilde{\varphi}_B) + a_E(\varphi_M, \widetilde{\varphi}_A), \end{align}\] showing that \[\begin{align} a_E(\varphi_M, \widetilde{\varphi}_B) = \frac{1}{2}. \end{align}\] On \(\Delta A'BC\) we have \(\widetilde{\varphi}_B + \widetilde{\varphi}_C + \widetilde{\varphi}_{A'} \equiv 1\), thus (implicitly) this holds on \(D_M^{(2)}\). So \[\begin{align} 0 &= & \int_{D_M^{(2)}} \nabla \varphi_M \cdot \nabla (\widetilde{\varphi}_B + \widetilde{\varphi}_C + \widetilde{\varphi}_{A'}) \stackrel{\eqref{eq:Dm2}}{=} \int_{D_M^{(2)}} \nabla \varphi_M \cdot \nabla \widetilde{\varphi}_C + \int_{D_M^{(2)}} \nabla \varphi_M \cdot \nabla \widetilde{\varphi}_{A'}\\ & \stackrel{\eqref{eq:MBC}}{=} & a_E(\varphi_M, \widetilde{\varphi}_C) + a_E(\varphi_M, \widetilde{\varphi}_{A'}), \end{align}\] showing that \[\begin{align} a_E(\varphi_M, \widetilde{\varphi}_C) = \frac{1}{2}, \end{align}\] as well. Therefore, \[\begin{align} && \widetilde{S}_{14} = \widetilde{{S}}_{41}= \widetilde{{S}}_{24} = \widetilde{{S}}_{42}=1/2, \end{align}\] which concludes the computation of the block \({{\mathbf{B}}}\) (recall that the other entries are \(0\)). All these results can also be obtained by classical trigonometric arguments. Furthermore, using the formulas 45 and 47 we obtain \[\begin{align} &&a_E(\varphi_M, \varphi_N) = a_E(\varphi_M, \varphi_P) = -\frac{\sin (2 \alpha)}{2 \sin^2\alpha}= -\frac{\cos \alpha }{ \sin\alpha} =-\frac{1}{\alpha} (1 + O(\alpha)). \end{align}\] The last entry needed is \[\begin{align} a_E(\varphi_M, \varphi_M) &=& \int_{{\mathcal{T}}_1\cup \dots \cup {\mathcal{T}}_5}|\nabla \varphi_M|^2 \\ &= & \frac{1}{\sin(\pi-2\alpha)} + \frac{\sin(\pi-2\alpha)}{2\sin^2 \alpha} + \overbrace{\int_{{\mathcal{T}}_4\cup {\mathcal{T}}_5}|\nabla \varphi_M|^2}^{O(1)}\\ &=& \frac{3}{2\alpha} (1+O(\alpha)), \end{align}\] which concludes the computation of \(\widetilde{{\mathbf{S}}}\). ◻