Remarks

Calculation speed increased 50% for calculating trial functions.

Automatically determine minimum polynomial order.
 Calculation speed increased >90% during calculation of residual integral. The way to do this is by implementing a new integration method. The new method and the old one show matched results.
Details
Dec 15, 2020 (Tuesday)

All boundary conditions implemented (clamped, simplysupported, free).

Checking whether results on edges satisfy known boundary conditions (unfortunately some of them do not).

Trying to debug the code. Starting from clamped cases and simplysupported cases.

One possible reason is the direction of the normal vector \(n\). The current code assumes the rotation of axis to be counterclockwise by default. This may cause the \(n\) vector to point opposed to the "inner area". This may be the reason why resulted slope on some edges have negative value. By drawing coordinates to visually / manually verify this theory. I believe this is the cause.

Trying to implement the correct \(n\) direction. Not easy. Will be working on this today the entire day I guess.

By the end of the day, the correct \(n\) direction has been implemented. Also for clamped and simplysupported edges, it has been verified by checking the corresponding moment / force / slope on edges. For both boundary conditions, visually the results on edges make sense. Free edge is still waiting for verification and possible debugging.
 More testing have been done on free edge. An interesting but disturbing result shows that for some combinations of clamped & free edges, the minimum order of polynomial gives the wrong / unreasonable result. It is known that increasing polynomial order will increase the accuracy, but this is unexpected since, by instinct, the minimum order of polynomial should also work somehow. Trying to fix this tomorrow.
Dec 16, 2020 (Wednesday)

Problem: For rectangle with boundary conditions as CCFC, the minimum solvable polynomial is with order of 10, which gives visually incorrect results. When polynomial order is increased to 11, still the visually incorrect results. When icreased to 12, the results start making sense. But the time it takes to calculate increased significantly.

To solve the problem, the first thing I want to do is to find a way to make the code run faster, especially when polynomial order is high. My reasoning for this is we probably need to use higher order polynomial to get more accurate results (visually correct ones) eventually. And it actually helps to make the code run faster anyways. So why not do it now.

One possible way to decrease the time it takes to run the code is to reduce the time it takes to do residual integral. This is tested with various cases that this step (residual integral) is the most timeconsuming step, as we increase the polynomial order.

A side note: function equationsToMatrix gives symbolic expressions (A, B), then using A\B+null(A) really takes a while to solve. Now A and B are converted to double first, then perform the calculation, which is significantly faster in some cases. A problem I had is A and B are symbolic, using null(A) gives me different results as using null(double(A)). After some research and tries, I found null(A)=null(double(A),'r'). The use of 'r' is explained here.

Also, trial function calculation is put in a loop, the initial polynomial order is small (N=5 for now). The loop checks if N will lead to nonempty trial function. If yes, we get the minimum polynomial order; if not, N is increased until we get nonempty trial function. This may be a temporary solution. It does some calculation until a solution is found, which takes time. But it's acceptable (for now).

Currently I am calculating residual integral by spliting the plates into slices of rectangles and calculate the area of each, then take the sum. I know there are other ways / builtin functions, but I have symbolic variables in the integral equations. I guess my way of calculating this is the easiest so far (I could be wrong, I will try different ways to do this integral, maybe it will be faster and more accurate). Anyways, spliting the plate into slices of rectangles, the more slices we have, the more accurate the result will be, but the longer it will take. I tried to use one case (rectangle with boundary conditions as CCFC, the minimum solvable polynomial is with order of 10) to calculate the time and integral results. I notice the followings: (1) the time it takes to calculate linearly increases with respect to the amount of slices; (2) the integral results (a coefficient) decrease exponentially with respect to the amount of slices (results closes to plateau after ~50 slices). I will verify this more with more cases, to determine what's the amount of slices we should use.
 The integral of residual function takes the longest time, as stated before. Within the calculation of residual integral, the part that calculate the integral expression (add all sliceintegral together) takes the longest time, about >90% of the calculation time. Now, its about to find another / faster way to solve this. I first tried to use different integral limits (ylimit), but it does not have any effect.
Dec 17, 2020 (Thursday)

First tried to write a new method to calculate residual integral. Hoped to be faster. It supposed to take out coefficient C and calculate the rest into matrix A and B, then use A\B to calculate C. Well, the method did somehow worked, but it took longer time to calculate. So I have given up on this method.

There must be a way to calculate double integral within a contour (ploygon), I am trying to find out how to solve this (I have done this by using small trapizoids, I want to find smarter / faster ways).

Trying another way to do residual integral. Using variable expressions for upper / lower limits. In this way, there will be a lot fewer slices (a lot fewer iterations of integral needed, save times), and more accurate (upper and lower bounds are now the exact boundary, not approximated rectangles). Results show much faster (80% improve).

Results validuated using rectangle with all edge simplysupported (existing analytical solution). Extremely good match, with error within 5%.
 Results also validated using other shapes and boundary conditions. Results on boundaries match the requested boundary conditions. However, results at other points can not be validated now, since there does not exist analytical solution. FEA is needed now, and this will be next step.
Dec 18, 2020 (Friday)

One problem in the code that I noticed this morning is that the new residual calculation does not make sense when polygon shape becomes random or the polynomial order increases. After degugging, I noticed that the way the code slice the polygon may give wrong solutions. The reason is the interation point of lines give much more accuracy (e.g. to the 50th decimal point) but visually they look the same. Therefore, the slices give wrong results for some random polygon. I tried to round up the intersection points to 10 decimal point and it works fine now. However, there must be another way to solve this problem (to neglect the high / redundent accuracy).

One thing I noticed during these trials: If the boundary conditions along all edges are the same, we only need minimum polynomial order to get visually correct results. However, if the boundary conditions are mixed, higher polynomial order may be necessary to get visually correct results. One more thing to deal with later.

Using the new residual integral calculation, the corfficient c is very close to the old residual integral calculation. (note: old residual integral calculation gives approximate c)
 Checked new residual integral calculation v.s. old way. The time they take is as follows:
Ploygon Shape  Boundary Conditions  Polynomial Order  Old Method Calculation Time (s)  New Method Calculation Time (s)  Speed Increase 

Right triangle  SSS  7(minimum)  17.7071  0.5290  97.01% 
Right triangle  SSS  8  104.3383  2.6930  97.42% 
Right triangle  SSS  9  434.9891  17.4255  95.99% 
Right triangle  SSS  10  1341.2622  94.3956  92.96% 
Random triangle  SSS  7(minimum)  9.7053  0.5142  94.7% 
Random triangle  SSS  8  47.6480  2.0723  95.65% 
Random triangle  SSS  9  309.8624  14.8110  95.22% 
Random triangle  SSS  10  1008.6408  51.5274  94.89% 
Random triangle  SCC  7(minimum)  16.5692  0.7590  95.42% 
Random triangle  SCC  8  96.7196  4.0388  95.82% 
Random triangle  SCC  9  413.7540  19.0235  95.4% 
Rectangle  SSSS  8(minimum)  11.0170  0.3224  97.07% 
Rectangle  SSSS  9  60.1881  1.0085  98.32% 
Rectangle  SSSS  10  276.4640  3.2389  98.83% 
Rectangle  SSSS  11  1032.9852  19.2665  98.13% 
Rectangle  CCCC  8(minimum)  10.1211  0.2797  97.24% 
Rectangle  CCCC  9  62.4365  0.8610  98.62% 
Rectangle  CCCC  10  286.7006  3.5423  98.76% 
Rectangle  SCSC  8 (minimum)  9.5832  0.3042  96.83% 
Rectangle  SCSC  9  62.0498  0.9841  98.41% 
Rectangle  SCSC  10  276.7912  3.4167  98.77% 
Rectangle  CSSC  8(minimum)  10.3034  0.2880  97.21% 
Rectangle  CSSC  9  61.3386  0.9279  98.49% 
Rectangle  CSSC  10  274.8270  3.4119  98.75% 
Random 4 edges  CSFS  11(minimum)  194.0299  10.4531  94.61% 
Random 4 edges  CSFS  12  1128.4976  67.9565  93.98% 
Random 4 edges  CSFS  13  3217.7059  217.5573  93.24% 
Random 5 edges  CCCCC  10(minimum)  18.3101  1.4104  92.3% 
Random 5 edges  SSSSS  13(minimum)  36.8810  3.2176  91.28% 
Random 5 edges  CSCSC  12(minimum)  28.7605  2.5393  91.17% 
Random 6 edges  CCCCCC  12(minimum)  26.0530  2.7281  89.53% 
Random 6 edges  SSSSSS  16(minimum)  1083.1359  132.9720  87.72% 
Random 6 edges  CSCSSC  15(minimum)  322.4051  38.8201  87.96% 
ToDo List

Free edge validation

Comparison of two plates.

Stress and strain calculation.
 The current implement works for convex polygons, and concave polygons with concaves vertically, but not horizontally.