I am currently studying this paper On Solving the IrregularlyShaped Plate Buckling Problem.

"Traditional finite element theory assumes a set of displacement functions that are nodally based and minimizes the potential energy with respcct to the nodal coordinates."

"The method of this paper provides a general displacement function with unknown coefficients and without regard to the element’s nodes."

"The only drawback to the present method is that it is limited to relatively simple plates defined by eight nodes and simply connected plates (no holes)."

Background knowledge 1: Gaussian quadrature. "Gauss quadrature allows integration of functions by simply evaluating functions at certain points, called Gauss points, rather than actually integrating using a technique such as Simpson's Rule or Romberg integration."

Background knowledge 2: Lagrange multiplier

Background knowledge 3: Subparametric mapping
 Background knowledge 4: Plane stress and plane strain
At the same time, I am researching other articles to solve what I don't understand in this paper.

"However, for many other boundary conditions, simple displacement functions like Equation (6.19) do not exist, and in some cases analytical, exact solutions analogous to Equations (6.21) and (6.31) have not been found. In those cases approximate solutions have been found using energy methods, which will be discussed in Chapters 8 and 9." See this article here: BUCKLING OF ISOTROPIC COLUMNS AND PLATES.pdf.
 RayleighRitz method for analysis (minimum potential energy theorm):
Precedures, Method Writeup

Step 1: Subparametric Mapping in FEA System. (I have no background in this, so it took a while to get the information I needed and some time to learn this knowledge.)

Goal: This is for transforming irregular shaped plane into square shape in another coordinate system. This will make integration easier.

See details in my note here.

I studied 'subparametric mapping'. From Cartesian domian into (\(\zeta\), \(\eta\)) domain.
 A figure of (\(\zeta\), \(\eta\)) domain, a quadrilateral element:

For a polygon, with number of edges no greater than 8, we can use 8node mapping method to change coordinates from Cartesian domian into (\(\zeta\), \(\eta\)) domain. I did not get the same result as that presented in the paper. I assume it is because the node sequence I use is different than the one used in the paper, which is not clearly stated. But my result also make sense. I will use the result from this paper in my later calculation.

The reason for doing this mapping is for easy integration in later steps.
 A good reference can be found here if you want to learn. Here is another one about quadrilateral element. Here is a useful website.
Step 2: Plane Stress, Strain Energy U Calculation.
 See details in my note here.
Step 3: Plane Stress, Virtual Work Energy T Calculation.
 See details in my note here.
Step 4: Solving Coefficients Using Minimum Potential Energy Method.
 Goal: get inplane displacement \(u\) and \(v\), in order to get plane stress and strain. This plane stress and strain will be later applied to calculate bucalking virtual work energy. The reason to do it this way is because, for irregular shaped polygon, it is hard to calculate loading conditions on edges.
 Detailed calculation steps / validations on Lagrange multiplier are shown below:
 See details in my note here.
Step 5: Buckling (Bending), Strain Energy U Calculation.

A full derivative of \(\frac{\partial \zeta}{\partial x}\), \(\frac{\partial \zeta}{\partial y}\), \(\frac{\partial \eta}{\partial x}\) and \(\frac{\partial \eta}{\partial y}\) can be solved using a method called implicit derivative and using Mathemetica, results can be found here.
 See details in my note here.
Step 6: Buckling (Bending), Virtual Work Energy T Calculation.
 See details in my note here.
Step 7: Solving Coefficients Using Minimum Potential Energy Method.
 Goal: calculate buckling load (coefficients).
 See details in my note here.
A Concrete Detailed Calculation
Summary:

{As of Jan 13, 2021} (Almost) All clear! Start working on implementation.

{As of Jan 11, 2021} There are not enough information in the paper that can lead me to solve the problem by writing down all necessary equations.

{As of Jan 11, 2021} The precedure seems to be limited to 8nodes only.

{As of Jan 11, 2021} The bending problem seems to be limited to rectangular only (if I understand correctly). And the expression of external load \(N_x\), \(N_y\) and \(N_{xy}\) seems to indicating the loading is quite specialized.
 {As of Jan 11, 2021} How to convert boundary conditions from one coordinate system to another is question to solve. In (\(\zeta\), \(\eta\)) domain, how to incorporate the BSc into the energy equation is still a question.