Element geometry in FEM for electrostatics In a study of mine, I am using tetrahedral elements to solve for magnetic flux density from a given current density distribution. I extract the mesh structure from COMSOL. Since it is much easier to implement, I have chosen tetrahedral elements. However, my advisor asked why I have chosen tetrahedral elements. How can I  justify this selection? I am using Gaussian quadrature to evaluate element-wise integrals. Thus, no help from that angle.

Tetrahedral elements are excellent for generating a mesh that captures polytopal domains exactly. Furthermore, tetrahdrons are affine equivalent, that is there is one reference tetrahedron $\hat K$, and all of the mesh tetrahedra $K$ can be obtained by applying an affine mapping $F_K:\hat K\to K$ to the reference element. Such maps are also crucial in the theory of polynomial approximation, and the affine equivalence is used to prove for example inverse estimates.
In contrast, one can use general polytopic meshes, but these are no longer affine equivalent in general. It is still possible to prove the well known results in polynomial approximation theory, but it is more work. Another thing to consider is that if you want to consider a continous piecewise polynomial space on a mesh consisting of arbitrary polytopes, you need higher order polynomials in order to maintain continuity (e.g. in 2D, the only continous piecewise affine function that is continous over a general mesh is in fact everywhere affine).
