Finite Element Method

The most general form of differential equation faced by PyNitride is

\[\begin{equation} \left[C^0(z)-iC^L(z)\partial_z-i\partial_zC^R(z)-\partial_z C^2(z) \partial_z \right] f(z) = w(z)b(z) \end{equation}\]

where the \(C\) are potentially matrices, \(C^L\) and \(C^R\) are transposes of one another. The \(C\) and \(w\) are material properties (ie they are well-defined on the mid mesh) whereas \(f(z)\) and \(b(z)\) are solution variables (well-defined on the node mesh). This can include \(b(z)=\lambda f(z)\) as an eigenvalue problem.

Let \(\Delta_i(z)\) be the linear interpolation of \(\delta_{z_i, z_j}\) on dependent variable \(z_j\). Choose the basis functions \(\Delta_i\) where \(i\), where \(i\) ranges over all the nodes except Dirichelet boundaries.

Now for each basis element, multiply the differential equation by that element and integrate. Treat the material properties as being step discontinuous to preserve sharp heterointerfaces. Wherever there is a material property to the right of a differential operator, use integration by parts to flip the differential operator onto the basis element instead. This gives us a \(N-d\) equations in \(N-d\) unknowns where \(N\) is the number of nodes and \(d\) is the number of Dirichelet boundary nodes. [*]

This equation can then be written \(Af=Mb\) or, for an eigenvalue problem, \(Af=\lambda Mf\), where \(A\) is the “stiffness matrix” and \(M\) is the “load matrix.” The stiffness matrix is given by

\[\begin{split}\begin{align} A_{i,i-1}= &\frac{C^0_{i-.5}dz_{i-.5}}{6}+\frac{i}{2}\left( C^L_{i-.5}+C^R_{i-.5} \right)-\frac{C^2_{i-.5}}{dz_{i-.5}}\\ A_{i,i+1}= &\frac{C^0_{i+.5}dz_{i+.5}}{6}-\frac{i}{2}\left( C^L_{i+.5}+C^R_{i+.5} \right)-\frac{C^2_{i+.5}}{dz_{i+.5}}\\ A_{i,i}= &\frac{C^0_{i-.5}dz_{i-.5}}{3}+\frac{C^0_{i+.5}dz_{i+.5}}{3}\\ &+\frac{i}{2}\left( -C^L_{i-.5}+C^L_{i+.5} + C^R_{i-.5}-C^R_{i+.5} \right)\\ &+\frac{C^2_{i-.5}}{dz_{i-.5}}+\frac{C^2_{i+.5}}{dz_{i+.5}} \end{align}\end{split}\]

And the load matrix can be written

\[\begin{split}\begin{align} M_{i,i-1}&= \frac{w_{i-.5}dz_{i-.5}}{6}\\ M_{i,i+1}&= \frac{w_{i+.5}dz_{i+.5}}{6}\\ M_{i,i}&= \frac{w_{i-.5}dz_{i-.5}}{3}+\frac{w_{i+.5}dz_{i+.5}}{3} \end{align}\end{split}\]

[*]

Note that, even at a Dirichelet boundary the next node in has support that extends all the way to the Dirichelet node, ie the edge of the solution domain, where it falls to zero. Thus the material parameters in that mid region at the edge will come into play even though that region is outside the range of node points actually used in deriving the equations below.

Contrast this to a Neumann boundary where the last node is at the edge of the solution domain, so there is no support of that basis element extending onto the outer side of the node. In the case of a Neumann boundary, the equations below which apply to each \(i\) have to be taken with the prescription that any term evaluated outside the node range is ignored. In the case of a Dirichelet boundary, that prescription does not matter since no terms will ever be outside the solution domain, because the solution domain is one node further than the set of equations.