where the matrix \(\mbox{$B$}{}^e\)  defines the gradient operator for a linear triangular element in particular. The details on the computation of this matrix can be found in the accompanying code. One of the advantages of this methodology is that the gradient is guaranteed to be tangent to the manifold, which is not the case for the traditional physics-informed neural network formulation [1]}. In a similar fashion, we can compute the Laplacian \(\Delta \hat{u}_i\)  at any given nodal location by predicting \(\hat{u}_j\)  at all the neighboring nodes, including the node of interest. Then, we can use the definition of the discrete Laplacian presented in equation REF  and approximate \(\Delta \hat{u}_i \approx \mbox{$A$}{}_{ij} \hat{u}_j\) .
