, where a is twice the average side length of the triangular finite element mesh (i.e., a = 2 • avg_side_length_poiss). This charge density is then stored in the variable rho_bottom_inside. While the numerical scheme for incorporating the charge density at the polar interface is relatively straightforward, implementing the surface charge for the semi-polar interfaces requires additional care since the Delaunay triangulation procedure does not create a symmetric grid. To this end, a rotation matrix is used in rho_semipolar_inside.m to first rotate the ( x, ỹ) coordinates of the finite element grid by 120 • counterclockwise so that one semi-polar interface now lies along the horizontal line ỹ = -c √ 3/6. A Gaussian having the same functional form as the one described previously (except with σ s replaced with -σ s /2) is then applied at the ỹ = -c √ 3/6 line within the domain -c/2 ≤ x ≤ c/2. The charge density is stored in the appropriate ( x, ỹ) locations within the variable rho_left_inside. The same procedure is repeated for the other semi-polar interface (except with a 120 • clockwise rotation) and the charge density is stored in the variable rho_right_inside. The three charge densities are finally added together and returned as an output variable by rho_semipolar_inside.m.